A/'S 6- 2d^ 


(NASA-CB- 15871 1) TIME. OPTIMAL CONTBOL OF A 
JET ENGINE USING A QUASI-BERHJTE 
INTERPOLATION MODEL M.S. Thesis {Notre Dame 
Univ.) 171 p HC A08/ME AOI CSCL 21E 


N79-25019 


Uaclas 


G3/07 23418 




ps 's\ 

, ft * ^ ^ Oil 

s' < ~x < s\ <k>ij 

&*. »s<b/ 


So 

%a e , % 


‘Department of 

ELECTRICAL ENgM^ERING 


UNIVERSITY OF NOTRE DAM E, NOTRE DAME, INDIANA 


TIME OPTIMAL CONTROL OF A JET ENGINE USING 
A QUA S I -HERMITE INTERPOLATION MODEL* 


John G. Comiskey 


May 1979 

Technical Report No. EE-791 


*This research was supported in part by the National Aeronautics and 
Space Administration under Grant NSG-3Q48. and also has been submitted 
in partial fulfillment of the requirements for the M.S.E.E. degree at 
the University of Notre Dame. 



ABSTRACT 


The scope of this work will be to make preliminary efforts to generate 
non-linear numerical models of a two-spooled turbofan jet engine, and sub- 
ject these models to a known method of generating global, non-linear, time 
optimal control. laws . The models will be derived numerically, directly from 
empirical data, as a first step in developing an automatic modelling proce- 
dure. 
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CHAPTER. I 


INTRODUCTION 

The scope of this work will be to make preliminary efforts to generate 
non-linear numerical models of a two— spooled turbofan jet engine, and sub- 
ject these models to a known method of generating global, non-linear, time 
optimal control laws. The models will be derived numerically, directly from 
empirical data, as a first step in developing an automatic modelling proce- 
dure. 

A hierarchy of models, including analytical and numerical models, will 
be established. The numerical models will be described in detail, and their 
step responses compared to those of the hypothetical jet engine from which 
they were derived. A method of generating time optimal control laws will be 
explained, programmed, and applied to the numerical models. Finally, these 
control laws will be tested, both on the models from which they were gener- 
ated, and on the hypothetical jet engine. 

This is the third in a series of similar works, whose ultimate goal 
is the development of an automated modelling method. Even though DYNGEN, 
an elaborate mathematical model, already exists, new models were developed 
for too reasons. First, DYNGEN uses too much cpu time to be called re- 
peatedly by an iterative method such as dynamic programming; a smaller, 
faster model is required. Second, DYNGEN assumes the role of a physical 
plant in this work, since access to a real jet engine is impossible. 

In. his paper. Basso [9] uses too methods of generating optimal control 
sequences. The first is the dynamic programming successive approximations 
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technique. This actually generates a control law, from which a control se- 
quence can be derived. The second is a modified Fletcher-Reeves conjugate 
gradient method. This method generates a control sequence that drives the 
system to the target in minimum time. The modification consists of the in- 
troduction of constraints into the original method. 

His findings were that both methods yielded similar results , and that 
the number of computations necessary to solve a problem increases geometric- 
ally with system order for dynamic programming, but only arithmetically for 
the conjugate gradient method. 

Longenbaker [1] applies the dynamic programming method to several models 
of the 1-100 engine. His models include several linear systems, and one 
non-linear, analytical system of differential equations derived from phy- 
sical and mathematical relationships among the state and control variables. 
Longenbaker concludes that the agreement between this analytical model and 
the DYNGEN simulator is not strong enough to justify gz’eat faith in the 
control law generated. 

In this paper, the same dynamic programming method is applied with a 
modification introduced to reduce cpu time, to several numerical, non- 
linear models of the F-100 jet engine. The conclusions are that better 
numerical agreement with DYNGEN is achieved by this numerical models than 
by Longenbaker ' s analytical model, with a much smaller expenditure of 
man hours. However, complex interpolation techniques cause these models 
to use extravagant amounts of cpu time. Either larger data bases and less 
interpolation, or a more economical technique like Basso's conjugate gra- 
dient method should be explored in future work. 



CHAPTER II 


TOO SPOOL TORBOFAN JET ENGINE MODELS 


2.1 Introduction 

A form was chosen for the system model, which isolates static and dy- 
namic portions of the system behavior, so that each of these can be modelled 
independently. Several methods of modelling these two portions were tried, 
resulting in a hierarchy of models. Each model was subjected to the same 
analysis for purposes of comparison. 

In this chapter, the system model form is derived, and the modelling 
methods outlined. These methods will be treated in detail in later chapters. 

2.2 Basic Approximation Approach [2] 

Now consider a method for obtaining nonlinear models. Let 

x = f(x,u) (2.2-1) 

with x an n vector and u an m vector denoting a dynamical system such 
as a jet engine, in which the state variables and parameters u remain pos- 
itive throughout the system operation and there is a function g(u) such 
that for each equilibrium point 

f(x,u) = 0 -<- ■> x = g (u) (2.2-2) 

The steady state system analysis involves the study of the function g(u). 

We propose to approximate the system (2.2-1) by 

x=A(x,u)[x - g(u)] (2.2-3) 

where A(x) is a square matrix which varies as a function of x. Notice 
that if Xp is an equilibrium point of (2.2-1), x^ = g(u D ), then a lin- 
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earization about this equilibrium point results in the linear system 

dx = ApSx + B^ySu (2.2-4) 

and a linearization of the approximating system (2) at x^ = g(u^) results in 

6x = A(x d )6x + [-ACx^) (u D )]6u (2.2-5) 

Hence, the linearization of (2) will match the linearization of (1) if and 
only if 

A( V ’ 5) • -5> If ( V ' b d <2 - 2 -« 

Also, if Ap is invertible, as is often the case for jet engine models, 
equation yields 

|a Gy » -y\ (2.2-7) 

These static and dynamic data are available from known algorithms [7], 
leaving only the choice of interpolation methods for generating non-linear 
models . 

2.3 Hierarchy of Models 

This work has resulted in the formation of a hierarchy of models, each 
a step in the development of an automated modelling method . They are clas- 
sified as follows : 

Model 0: The actual F-100 type engine (hypothetical) 

Model 1: The DYNGEN [6] simulation program, coded with data presumed' to 

have been taken from experimental measurements on Model 0. This 
model solves 16 nonlinear differential equations and uses data 
maps and thermodynamic tables which cannot be expressed analy- 
tically . 
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Model 2: The Longenbaker [1] model, a 5th order, nonlinear, analytical mod- 

el. It includes the 5 state differential equations which govern 
the dynamical behavior of the system, along with 20 algebraic 
equations which express the relationship between various engine 
variables. This model is discussed in detail in [1]. 

Model 3A: The linear affine power law model, which is a fit of steady state 
data to a selected form with linear, nonlinear and constant terms. 

Model 3B: The straight linear affine model, generated in the same manner 
as 3A, without the non-linear terms, to serve as a comparison. 

Model 4: The Quasi-Hermite interpolation model. Also a fit to steady 

state data, this model employs value and derivative matching 
over a two dimensional subset of the state space. 

Models 3 and 4 will be outlined briefly here, and detailed in later 
chapters. 

2.4 Linear Affine Power law Model [2] 

This model approximates the system by interpolating values of A(x,u) 

from values of the matrix at two data points, and by generating values, for 
g(u) by a fit of the form: 

C 3i C 4i 

g i (n) = c li u 1 + c 2i u 2 + c 5i u 2 u 2 + c &i , i = 1, . . . ,5 (2.4-1) 

to the same two data points. 

2.5 Quasi-Hermite Interpolation Model 

This model approximates the system by interpolating values of A(x,u) 
from values of the matrix at five data points, and by interpolating values 
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of g(u) from 15 data points using a two dimensional adaptation of Hermite 
interpolation. This model represents new work for this thesis. 

2 . 6 Summary 

Having chosen to express the model in the form: 

x = A(x,u) [x - g(u)] (2.2-3) 

it remains to derive the function g(u) and the matrix A(x,u) to corres- 
pond to empirical data. The function g(u) represents a mapping from the 
control space U into the state space X which yields steady state values 
for given controls. Empirical data available (i.e. DYGABCD output) includes 
both steady state values, and derivatives at those points with respect to 
the various inputs. It is desirable to choose functions which match as 
many of the available data as possible. 



CHAPTER III 


LINEAR AFFINE POWER LAW MODEL 3 


3.1 Formation of A(x,u) [2] 

Values of A(x,u) are interpolated from the values of Ap = ACx^jU^) 
and A^ = A(x^,u^ ) in the following manner: 


x. - 


A(x,u) = A w diag (~~~^ „ J 0 + diag(^ — 

Dj Wj D_. Wj 

where diag(-) is a diagonal matrix which causes the j ^ column of 
be interpolated linearly between the jth columns of A^ and A^ 
as the interpolation variable. 


(3.1-1) 

A(x) to 

with x . 

1 


3.2 Approximation of g(u) [2] 

The parameter vector u is presumed to be made up of physical control 
variables, and parameters such as fuel flow and nozzle area. The equilibrium 
function is to be approximated in a manner such that both the equilibrium 
values and the linearizations of the approximating system (3) match those 
of system (1) at both x^ and x^. This requires then that 


g(u D ) = Xp gCu^) = x w 


(3.2-1) 


and also 


f <v ■ V b d - 15 ( v - -\\ 


(3.2-2) 


The method we propose here is to approximate each scalar component g^(u) 

of g(u) by a linear affine power law form 

c c c 

\ , , . m+1 m+2 m+m , m n r>\ 

C 1 U 1 + ’” C m U m + C 2mFl U l U 2 * * * U m + C 2mP2 (3.2-3) 


for which the jth partial derivative is 
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£fi 

6u. 

3 


C. + C ')_L1 C ^L' 
j 2m+l m+j 


c c 

m+1 m+m 
u. . . . u 
1 m 

u. 


(3.2-4) 


Now, if the variables are normalized and scaled such that 


^ — (1,1, «.*,!) ^ (a, a, « • . ,a) a 

then, the conditions of (11) and (12) can be put in the form 


(3.2-5) 


6g i 

k j " 6Uj “ C j + C 2nrf-l C m+j 

5g i ^Tn+r 1 

Vj = 6^7 ® = C j + C 2mbl C m+j a 


v 2m+l 


= %tt) - Ifi. + + C; 


k 2^+2 * S i (a) ’ a ^i + 


2m+2 


*Vj . 

p Q X p 

j ' 2m+l 2nt+2 


(3.2-6) 


and summing the first two of these over j yields 


Z k J ’ + c 2mfl Krf-j 

£c .~1 

u ^1+ j 


^Sn+j ^ C j + C 2m+l a 
k 


2m+l 


= ^ C j + C 2m+1 + C 2m+2 


h-. 

L m+j 


k 2m+2 = a ^ C j + C 2mfl a “*' J + C 


2m+2 


(3.2-7) 


which is of the form 


S 1 " r l + r 3 r 2 


s 2 = r l + r 3 r 2 a 


V 1 


s 3 " r l + r 3 + r 4 


s, = ar, + r~a + r. 

4 13 4 


(3.2-8) 
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which, incidentally is the m=l condition also. This set of transcendental 
equations is solved numerically for r^, r^, and (3.2-6) is then 

used to solve for each c^. . In the event that (3.2-8) has no solution, a 
best fit is made on the second equation by varying while the other con- 
ditions are satisfied exactly. 


3.3 Computational Algorithm [2] 

In this section, we present an algorithm which serves to automate the 
process of finding a nonlinear model for a system 


x = f(x,u) 


(3.3-1) 


to be approximated from x^u^x^u^, by a normalized system. 
Hie algorithm will automatically perform the normalization and, hence, ac- 
tually approximate the system 


x = f(x,u) 


(3.3-2) 


where x. 

x 


where 


x. 


:^/xp , u. = u./ix^ . The approximating system is of the form 

1 i J J 3 


x = A(x) [x-g(u) ] 


V “*i „ *i “V 

A<*> - h dia8 + hi dia8 


x x 


i x 


and 


h 


“ I 


where 


c%* + ci _ TT u* m 3 4- 

“33 2 m+l ■ 3 2 mf 2 

3 3 


u* = a ,<L + 3 . . 
I 3 j 3 


Algorithm I . 

1. Input: VVVV m ’ n ’ a ’ £, VVV B W 

where m = number of controls 


(3.3-3) 


(3.3-4) 


(3.3-5) 


(3.3-6) 


n = number of states 
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2. Calculate: 


Ajj = diag(l/x D ) /^diagCXp ) 

X X 

\ = diagCl/Xjj ) A^diagCXp ) 
i i 

A. 

B d = diag(l/x D ) B D diag(u D ) 
x i 

A 

B w “ diagCl/Xp ) B w diag(u D ) 


Calculate: 

« - (l-a)^ /(u -^ > 

111 

S " (au D.-"w. )/(u D.-“w. ) 

3 3 3 3 

Calculate: 


3 ~ 1 » • ♦ • 


- <-*» Vtj 


k = 1 = 
2m+l 


*D. 


m+j 

Calculate 


*WH " KViJ k Lrf-2 “ V /jt D. = V 


3 1 j • • • jIQ 

x X) * • • jH. 


m 


1 3=1 


S 1 = k 1 
2 2m+l 


m 


s 3 " = T k 1 
2 


x . x 
S 4 k 2m+2 


X * 1 y 1 ■ • j XI 


6. Go to Algorithm II. 


Send: 


^1* S 2 > s 3’ s 4’ ^ 5 ^ 
Receive: r^, r^, r 3 , r^, y 


x lj . « * jii 


7. Calculate: 
i 

C 2m+1 


= r„ 


"2m+2 


j = 1, . . . ,m 
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k* . - k* + y . 

c! = k" 

1 v 1 3 3 

r^Ca -1) 


i i 
r o c . • 

3 uri-j 


8 . Output : 


c l’"*’ c 2m+2 


a. , P. 
i 3 


X D.> V 

i l 


V \ 


Algorithm II. 


1 * Inpu t. 3^ 


2. Calculate: 


S 3" S 4 

P 1 = S ! " "3T" 


s 3" s 4 

P 2 S 2 “ 1-a 


Minimize by line search: 


x-1 a -1 


P 2 " P 1 


for -10 < x < 10, x ^ 0, x # 1 


4. Calculate: 


2 r 3 r 

a -1 
r 2 " a-1 

r 2 

s - s, + r,(a -a) , r - 1 ' 

r l = 1=1 Y = m (S 1" S 2 + r 2 r 3 (a “ 1)} 
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r 2 \ 

- as^ - r^Ca -a) 

r 2 = 1 ^ 

5. Return to Algorithm 1.6 

3.4 Straight Linear Affine Model 

As a check that the Power law term has significant effect. on the function 
g(u), a straight linear affine approximation to g(u) was generated. This 
model is then subjected to the same analysis as models 3A and 4. 


3.5 Numerical Results 


[ 2 ] 


The algorithm of the previous section was applied to data obtained 


using DYNGEN with x^ and specified as in Section 2. An off-design 


point was obtained using u^ = (.72727, .72727), with the resulting norm- 


alized state x^ = (.9000, .7897, .7381, .9401, .9454). The normalized A 
and B matrices are 






-3.8 

-1.277 

2.067 

-1.152 

2.748 

-5.39 

1.585 

-1.991 

377.9 

49.51 

-264.9 

86.807 

31.26 

139.3.9 

-6.269 

-88.69 

-176.5 

23.91 

-10.27 

-37.4 


-4.744 

-1.3888 

3.2468 

-1.4591 

.82.86 

-26.726 

2.5585 

-1.8609 

475.73 

137.55 

-328.91 

27.791 

-50.103 

110.91 

63.188 

-116.69 

-186.77 

-67.682 

-41.681 

24.586 


1.448 


-.00259 

.3553 

1.071 

A 

.2116 

-.31618 

78.91 

B w = 

12.54- 

-13.774 

27.83 

-.6201 

-99.3 

246.7 


157.78 

6.84 






(3.5-1) 


1.1969 


-.04546 

.0013 

.45548 

A 

.0086 

-.0121 

91.495 


2.434 

-.613 

8.2883 

.67865 

-97.467 

243.23 


203.44 

.64755 



U— 

— 


(3.5-2) 


Using the parameter value a = .7, the c^ coefficients which specify the 
equilibrium functions §(u)^ as in Section 3 are given by the matrix . 
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.24267 

-.00218 

1.90082 

8.09916 

.02864 

.73088 



1.01593 

.85407 

.89872 

.66919 - 

.81879 

-.05121 


c = 

.73445 

.10133 

6.90586 

3.09409 

.011495 

.15272 

(3. 


.77234 

-.35905 

2.45867 

2.87415 - 

075198 

.66191 



.39503 

-.27262 

-3.44682 

13.4468 

.01838 

.85921 

' «*■ 


This matrix together with 1 

the values 

a=l.l and 

8=0.1 

and the 

matrices 


Ap and completely specify Model 3A. 


Another model which we will call Model 3B is easily obtained by using 


a linear 

affine approximation 

to j 

!(u) 

such 

that 

£ ( V = V fiC V = V 

Model 3B 

is specified by a = 

-1 
e > 

a = 

2.31778, 8 

= -1.31778 and the co~ 

efficient 

matrix 









.1553 

.0028 

1.0 

1.0 

0. 

.8418 




.1619 

.1707 

1.0 

1.0 

0. 

.6674 



G = 

.5351 

-.1208 

1.0 

1.0 

0. 

.5857 

(3.5-4) 



.5878 

-.49313 

1.0 

1.0 

0. 

.9053 




.2962 

-.2099 

1.0 

1.0 

0. 

.9137 




CHAPTER IV 


QUAS I-HESMITE INTERPOLATION MODEL 4 


4.1 Introduction 

A known method for matching a polynomial to the values and derivatives 
of a function at several points is Hermite interpolation. However, this 
method is formulated in general only for the one dimensional case. [3] Some 
works exist which apply this method to an n dimensional case, hut only 
under certain narrow restrictions. [5] Ihe single variable case, .its -restricted 
application in n dimensions , and a modified application in two dimensions , 
are discussed in this chapter. 

4 . 2 Hermite Interpolation for a Single Variable 

This presentation of the Hermite interpolation method is drawn from 
Hildebrand. [3] His notation is preserved, as closely as possible, here 
and in the resultant computer program. 

12 Til 

If the values of g(u) are known at m points, u = u , u ,...,u , 
define: 

fr (u) = (u - u 1 ) (u - u 2 ) . . . (u - u™) (4.2.1) 

and: 

£ i (u) = .. (4.2.2) 

(u-u ) ir"(u ) 

with the properties : 

tt(u^) =0 j = l,...,m (4.2.3) 

and : 

= S I = l,...,m j = 1, . . . ,m (4.2.4) 
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where 6.. is the Kronecker delta (5.. =1. 6.. =0 for all i ^ i). 

2.3 ii ij J 

With these defined, the polynomial of degree m-1 which takes on the 

values gCu 1 ) , g(u 2 ) » . . . , g(u m ) can be expressed as: 

y(u) = l £ k (u)g(u k ) (4.2.5) 

k=l 

12 in 

Suppose both g(u) and g^u) are known for u = u , u , ...,u , it 
is possible to determine a polynomial of degree 2m-l with these values 
and derivatives. We shall assume this polynomial is expressible in the 
form: 

y(u) = £ h k (u) g(u k ) + l h k (u) g"(u k ) (4.2.6) 

k=l k=l 

where h X (u) and h^u), i = l,...,m are polynomials of maximum degree 
2m-l. The requirement y^) = g^) will be satisfied if: 

h 1 (u^ ) =6.. and h X (u j ) = 0 (4.2.7) 

3-3 

and the requirement y'(u J ) = g'(u~*) will be satisfied if: 

h' X (u^) - 0 and h" 1 ^) - 6.. (4.2.8) 

Since £ (u) is a polynomial of degree m-1 which satisfied (4.2.4), 
then [^(u)] 2 is a polynomial of degree 2m~2 which satisfies (4.2.4) 
and whose derivative is zero at u/^ when i ^ j. So if h X (u) and h X (u) 
are polynomials of degree 2m-l, then: 

^(u) = /(uH-e^u)] 2 and h^u) = s^iOttV)] 2 (4.2.9) 

where r x (u) and s x (u) are linear functions of u, so that (4.2.7) and 
(4.2.8) will be satisfied when i ^ j . These four conditions, when i = j, 
then yield: 
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/(u 1 ) = 1 r^Cu 1 ) +.2 -^(u*) = 0 
s^Cu 1 ) = 0 s^Ci/) = 1 


(4.2.10) 

(4.2.11) 


from which follows: 

r^(u) = 1-2 £" i (u i ) (u - u 1 ) and s^(u) = (u - u^) (4.2.12) 

So, by combining (4.2.6), (4.2.9), and (4.2.12) we obtain the desired poly- 
nomial in the form: 


where : 


and: 



(4.2.13) 


(4.2.14) 


r*(u) = 1-2 £^ i (n i )(u - a*) and s i (u) = (u - u*) . (4.2.15) 

This result is known as Hermite T s interpolation formula, or the formula for 
osculating interpolation. 


4.3 Problems of Hermite Interpolation in ri Dimensions 

If Hermite Interpolation were to be applied in n dimensions, the 
task would be to determine m sets of n+1 polynomials with properties sim- 
ilar to h and h. Specifically, if we assume that the desired polynomial 
can be expressed in the form: 

y(u) = l h (u) g(u ) + l l h jk (u) (u k ) (4.3.1) 

k=l j=l k=l du J 

then these polynomials must have the properties: 

h 1 (u^ ) = 6 ± j and h lk (u^) =0 (4.3.2) 


and: 
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Sh 1 , k. _ , 6h jk , k. . 

(u ) = 0 and (u ) = 


(4.3.3) 


6u 


Su 


3 1 

corresponding to the conditions (4.2.7) and (4.2.8). However, the further 
condition: 


-i£ 

— — (u k ) = 0 for all i ^ k 

du. 

2 


(4.3.4) 


must also be satisfied. This final condition cannot be satisfied by the 
polynomials described in the previous section. In [5], Niijima treats a 
special case, in which the existence of certain orthogonal polynomials al- 
lows the application of Hermite interpolation to carefully chosen data in 
two dimensions. However, this method is not universally applicable to ar- 
bitrary data. 


4.4 Quasi-Hermite Approach for Two Controls 

Given that no general method of Hermite interpolation in taro variables 
was found, the following adaptation of the one dimensional case was applied. 
The value of control u^ (nozzle area) was held constant at the design 
point value, and Hermite interpolation was applied to a set of points gen- 
erated by varying u^ (fuel flow) . Both values , and derivatives with re- 
spect to u^ were matched at these data points. Then, for each value of 
u^, a value A was chosen, and control u^ was varied by this amount, 
both plus and minus. A function was then chosen to match values at these 
new points, without altering the function at the original points. The re- 
sulting polynomial is of the form: 

y(u) = \ hk ( u i) ek ( u 2 ) + l 

k=l k=l 


• (4.4.1) 
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where: 


z k 1 , .k. f k 1 ,k, 

k g(u l» u 2 + A ' ” 8C«-i » “ A ' 

8 (u 2 J = 1 + [ 73 71737 u 2 

2 A g(u lS u 2 ) 


, k 1 .k. , k 1 • .k. „ , k 1. 

g(«l» u 2 + A ) + g(u^,u 2 - A ) - 2g(u 1} u 2 ; 2 

+ [ : “ro"-;"— i- — <u 2 ) 

2(A ) g(« u ) 


This function has the property that: 


k 1 

0 (u 2 ) = 1 when u 2 = u 2 


0 1<: (u 2 ) = g(u 2 ,u 2 + A k )/g(u k ) when u 2 = u 2 + ^ 
0 k (u 2 ) = g(u k ,u 2 - A k )/g(u k ) when u 2 = u 2 “ 


and since h (u ) = 6.., the resultant polynomial will match values and de- 

■^-3 

rivatives with respect to u^ along the u 2 = u 2 line, and values at u 2 = 

1 , A k 
u 2 ± A . 


4.5 Formation of A(x,u) 

Having chosen an approximation to g(u), it remains to choose a method 
for interpolating values for A(x,u) to complete the model. Lagrangian 
interpolation was used to match values only at three points along the u 2 = 
u 2 line, and at u = (u^, u P + A 3 ). The results of this approach are em- 
bodied in the following equations. First define: 


A1 

= A(x,u) 

at 

u = 

1 

u — 

Cuj, 

V’ 

x = 

gCu 1 ) 

A3 

= A(x,u) 

at 

u — 

3 

u = 

(u l> 

3 * 

u 2 ). 

x = 

g(u 3 ) 

A5 

= A(x,u) 

at 

u = 

5 

u = 

(u l* 

5. 

u 2 )> 

X = 

g(u 5 ) 

AP 

= A(x,u) 

at 

u = 

P 

u = 

(u r 

u 2 + 

A 1 ) ,x = 

g(u P ) 

AM 

= A’(x,u) 

at 

u = 

M 

u - 

Cu r 

1 

U 2 ~ 

A^) ,x = 

g(u M ) 


(4.5.1) 
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and: 


then: 


FR1 = [ (x 1 - x^) ^ - x^) ] / [ (x* - x^) (x^ - x^) ] 

FR3 = [ (x 1 - xj) (x 1 - x^) ] / [ (x^ - x*) (x^ - x^) ] 

FR5 = [(x x - x i)( x 1 ~ “ X J)( X 2 “ x i^ 

FKP = [(u 2 - u 2 ) (u 2 - (uj - A 1 )) ]/[2(A 1 ) 2 ] 

FR = t(n 2 - (u 2 + A 1 )) (u 2 - (u 2 - A 1 ) ]/ [(-(A 1 ) 2 ] 
FRM= [(u 2 - (u 2 + A))(u 2 - u^l/UCA 1 ) 2 ] 

A ±j (x,u) = [FR1 (Al^j) + FR3(A3 ±j ) + FR5(A5 ± ^) ] 

x [FEP(AP../A1..) + FR + FRM(AM. ./Al. .)] 
ij id ij 13 


( 4 . 5 . 2 ) 


6 . 5 . 3 ) 



CHAPTER V 


MODEL RESPONSE COMPARISONS 


5.1 In troduc tion 

Before subjecting the models to the Dynamic Programming Algorithm, 
some effort was made to examine their closeness of fit to DYNGEN data. 

Steady state values of models 1 and 4 are compared, and fuel flow step re- 
sponses of models 3A, 3B, and 4 are plotted in comparison to DYNGEN responses. 
A description of the step response program is also included. 

5.2 Steady State Comparison of Models 1 and 4 

The function g(u) represents a mapping from the control space into 
the state space, relating fixed controls to steady states. It is not only 
useful in the model form: 

x = A(x,u) [x - g(u)] (5.2-1) 

but should also approximate the operating line of the plant. 

Such a comparison is made here between g(u) for model 4 and the 
DYNGEN simulator. Nozzle area was held constant, as fuel flow was varied 
from 9.0 to 1.1 by 0.02. All values are normalized. Percentage error is 
also computed, and shows the model’s excellent agreement in its range of 
accuracy, and rapid deterioration outside that range. 


20 




21 


Table 5.2-1 
X(l) = NC 


fuel flow 


DYKGEN 


Model 4 


% error 


0.90 

.97275 

.97288 

0.01 

0.92 

.97761 

.97790 

0.03 

0.94 

.98326 

.98326 

0.00 

0.96 

.98887 

.98892 

0.01 

0.98 

.99445 

.99461 

0.02 

1.00 

1.0000 

1.0000 

0.00 

1.02 

1.0051 

1.0051 

0.00 

1.04 

1.0102 

1.0113 

0.11 

1.06 

1.0152 

1.0230 

0.77 

1.08 

1.0201 

1.0513 

3.06 

1.10 

1.0244 

1.1195 

9.28 
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Table 5.2-2 


X(2) = NF 


el flow 

DYNGEN 

0.90 

.97132 

0.92 

.97883 

0.94 

.98427 

0.96 

.98961 

1.00 

1.0000 

1.02 

1.0046 

1.04 

1.0091 

1.-06 

1.0136 . 

1.08 

1.0180 

1.10 

1.0221 


Model 4 

% error 

.97099 

-0.03 

.97817 

-0.07 

.98425 

0.00 

.98948 

-0.01 

1.0000 

0.00 

1.0011 

-0.35 

.98046 

-2.84 

.89094 

-12 . 10 

.62787 

-38.32 

-.013230 

-101.29 
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Table 5.2-3 
X(3) = P4 


fuel flow 


DYKGEN 


Model 4 


% error 


0.90 

.92038 

.92042 

0.00 

0.92 

.93623 

.93633 

0.01 

0.94 

.95225 

.95225 

0.00 

0.96 

.96821 

.96824 

0.00 

0.98 

.98413 

.98434 

0.02 

1.00 

1.0000 

1.0000 

0.00 

1.02 

1.0148 

1.0125 

-0.23 

1.04 

1.0295 

1.0136 

-1.54 

1.06 

1.0441 

.98374 

-5.78 

1.08 

1.0587 

.88135 

-16.75 

1.10 

1.0727 

.62822 

-41.44 
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Table 5.2-4 


X(4) = P7 


lel flow 

DYNGEN 

0.90 

.94118 

0.92 

.95358 

0.94 

.96532 

0.96 

.97697 

•0.98 

.98853 

1.00 

1.0000 

1.02 

1.0107 

1.04 

1.0213 

1.06 

1.0319 

1.08 

1.0425 

1.10 

1.0527 


Model 4 

% error 

.94109 

-0.01 

.95340 

-0.02 

.96531 

0.00 

.97693 

0.00 

.98856 

0.00 

1.0000 

0.00 

1,0081 

-0.26 

1.0013 

-1.96 

.94911 

-8.02 

.78384 

-24.81 

.37203 

—64 . 66 
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Table 5.2-5 


X(5) = U4 


iel flow 

DYNGEN 

0.90 

.96625 

0.92 

.97304 

0.94 

.97992 

0.96 

.98670 

0.98 

.99339 

1.00 

1.0000 

1.02 

1.0074 

1.04 

1.0147 

1.06 

1.0219 

1.08 

1.0290 

1.10 

1.0365 


Model 4 

% error 

.96624 

0.00 

.97308 

0.00 

.97992 

0.00 

.98665 

-0.01 

.99320 

-0.02 

1.0000 

0.00 

1.0089 

0.15 

1.0248 

1.00 

1.0573 

3.46 

1.1231 

9.14 

1.2494 

20.54 
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5.3 Program Layout 

The method chosen for generating time response data was a Euler in- 
egration with a user varied time step. After specifying initial controls, 
the user provides a control sequence of time step, duration (in iterations), 
and controls. This structure allows the user to provide smaller time in- 
crements for the steeper portions of the response, and to alter the con- 
trols during the response. The step response program creates a file of 
time-state n- tuples, which are plotted against similar DYNGEN data by an- 
other program. 

5.4 Euel Step Response for Models 3A, 3B, and 4 

Each of the three models was subjected to a fuel flow step from 0.8 
to 1.0, and the response plotted against the same response by DYNGEN. 

These graphs show that all three models match DYNGEN closely, but that 
model 4 is a better fit than either 3A or 3B. 
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■Figure 5.4-1 





STEP RESPONSE OF MODEL 30 VS. QYNGEN 
FOR FUEL FLOW STEP FROM 0.8 TO 1.0 
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Figure 5.4-4 
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Figure 5.4-7 
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Figure 5.4-10 










CHAPTER VI 


SUCCESSIVE APPROXIMATION DYNAMIC PROGRAMMING 

6.1 Introduction 

The successive approximation dynamic programming method is described in 
detail by Longenbaker [1], This chapter will state the problem, and describe 
refinements made to the previous software. In addition, the general struc- 
ture of the software will be discussed, and repitition of an example from 
Longenbaker will serve to verify its accuracy. 

6.2 Time Optimal Control Problem [1] 

The necessary first step is to reformulate the models as discrete time 

systems. Let 

x(t + At) = x ( t) + At-f(x(t), u (t) ) (6.2-1) 

represent the system with starting time k and terminal time N. It is 
understood that 

f(x(t), u(t)) = Ax(t) + B u(t) (6.2-2) 

for linear models. Let x(k) be the starting state and let the terminal 
time N be defined as the first instant at which the system state reaches 
the designated target set S. All x(t) are eX, the state set. The per- 
formance index 

N-At 

j(x,u) = l At (6.2-3) 

t=k 

t = k, k + At, . . . . ,N-At 

is to be minimized with u(t) e U, the control set, and _u defined as the 
control sequence: 
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u_=u(k), u(k + At) , . . .u(N-At) (6.2-4) 

Furthermore, the minimization is subject to hard constraints of the form 

y i (x(t), u (t) ) £ c.. (6.2-5) 

6.3 Refinement of the Dynamic Programming Method 

The dynamic programming technique requires that the initial estimates 
of the cost function be greater than or equal to the final values. Further- 
more, the method guarantees that the output of each iteration meets the same 
requirement. There are, however, two sources of error in the intermediate 
outputs: error due to a non-optimal choice of controls, and error due to 

overly large initial estimates. The former can only be eliminated by searches 
of the entire control space, but the latter can be eliminated by repeated 
iterations on a fixed control law. 

Each control law, whether optimal or non-optimal, yields a set of cost 
values, which can be determined by the same successive approximations tech- 
nique. By the definition of optimal, the cost values generated by any con- 
trol law are guaranteed to be greater than or equal to the corresponding 
optimal values, and therefore qualify as input to an iteration of the gen- 
eral dynamic programming technique. 

To effect this change, a subiteration loop repeats the successive ap- 
proximations technique on a fixed control law until the values converge on 
the actual costs for that fixed set of controls. These values are then 
used as input to another iteration, including a search of the entire con- 
trol space for a better control law. In this manner, a large number of 
control space searches are eliminated and a similarly large amount of c.p.u. 
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time is saved. 

6.4 G.eneral Program Structure 

Since the dynamic programming method was to be applied to several dis- 
similar models, an effort was made to isolate the two portions of the sys- 
tem, the dynamic programming method and the particular model, as much as 
possible. To that end, all the models were constructed as sets of subrou- 
tines, with matching subprogram names. The main program then refers to these 
routines to access model dependent quantities . 

A brief description of the main program routine follows : 

MAIN: The dynamic programming algorithm: reads data describing 

the state and control spaces, initial cost estimates, and 
iteration control values; performs iterations and subiter- 
ations until convergence occurs or iteration counts are ex- 
hausted; records results on disk and printer. 

SPIRAL: This subroutine computes the indices of the next 

statespace point, spiralling out from the target, as de- 
scribed by Longenbaker [1] . 

V: This function subroutine interpolates new values of 

the cost function from current state and cost data. 

NEXTX: This subroutine computes x(t 4- At) and tests 

the new value against state space and output constraints . 

The routines refered to by the dynamic programming algorithm, and there- 
fore required in all models are: 

INIT: This subroutine performs any initialization nec- 
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essary for operation of the model. 

XDOT: This subroutine computes derivative data as a 

function of state and control. 

OUTPUT: This subroutine computes output data as a function 

of state and control, which may then be used in constraint 
tes ting . 

COMPLT: This subroutine computes values for states not 

specified by the dynamic programming algorithm, should 
the order of the model be greater than two. 

With this structure, all that need be done to change models is to con- 
catenate the FORTRAN source for the dynamic programming algorithm and the 
desired model, and provide reference to any data sets the model may require. 
For more information concerning the operation of the algorithm, refer to 
the commented listing in appendix F. 

6.5 Verification of Longenbaker 1L2 Unconstrained 

Since new dynamic programming software was developed to increase gen- 
erality and operating efficiency, it is necessary to verify proper perfor- 
mance. To this end, comparison was made between the results of the current 
software and that of Longenbaker for a linear example. In the unconstrained 
case, cost values differ by no more than 0.0001 seconds, which can be con- 
sidered insignificant. The control laws differ in only a few cases, and 
then only by a single control space increment. Given the inherent inaccuracy 
of the control space quantization, this too is insignificant. 
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FIGURE 6.5-1-B. Model 1L2 (Uncon.stra.ined) - Cost as per Longenbaker [1] 
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6. 6 Verification of Longenbaker 1L2 Constrained 

There is considerable disagreement between the results of the two pro- 
grams, giving rise to the question:' Which, if either, of the results are 
correct? Without extensive tracing of program operation, this question can- 
not be completely resolved, however, there is some evidence to motivate a 
choice of the new software. 

Consider the case of x(l) — 1.02 and x(2) — 1.00. Both programs 
select identical optimal controls, u(l) = 0.50 and u(2) = 1.20. The same 
controls applied to a state adjacent (in terms of the state space quantiza- 
tion) to the target should yield the same cost value. This is not the case, 
with the Longenbaker software giving 0.0735 seconds and the new software 
giving 0.0606 seconds. 

Note, however, that the unconstrained case yields precisely the same 
control (u(l) = 0.50, u(2) = 1.20), and both the Longenbaker and new soft- 
ware give cost values of 0.0605 seconds. With or withcut constraints, the 
same control applied to the same state adjacent to the target should yield 
the same cost. With three of four calculations in agreement, it can only 
be concluded that the fourth is incorrect. The source of this disagreement, 
especially in light of the agreement with the unconstrained case, is func- 


tionally indeterminable. 
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Model 1L2 (Constrained) Optimal Control Lav as per Longenbaker [1]. 
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FIGURE 6.6.1-B. Model 1L2 (Constrained) - Cost as per Longenbaker [1] . 
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CHAPTER VII 


OPTIMAL CONTROL LAWS 


7.1 Introduction 

Before subjecting each model to the dynamic programming algorithm, care 
must be taken in choosing state and control space parameters, to insure the 
accuracy of the results. The target point and state space quantization must 
be chosen so that the entire state space window lies within the model's range 
of accuracy. The steady state controls for the target must also be known, 
since they will be needed during simulation. The third of model 4's five 
data points (u(l) = 0.'8, u(2) = 1.0; x(l) = 0.948434, x(2) = 0.928751) was 
chosen as the target and 0.01 as the quantization, to satisfy these require- 
ments . 

Control space limits must also be chosen to guarantee accuracy. Only 
controls within the range of those used to generate the data points can be 
used with certainty. For models 3A and 3B, 0.74 ^ u(l) <.1.0 and 0.74 £ 
u(2) 1.0 were used; for model 4 the control space limits are 0.74 _< u(l) 

j< 1.0 and 0.87 u(2) j< 1.13. These limits correspond closely to the con- 

trols used in data point generation.* 

7 . 2 Model 3A Constrained 

The extremal controls appear almost exclusively throughout three quar- 
ters of the state space windoxj, indicating that system performance is re- 
stricted by the control space limits. These limits may not be relaxed, 
however, since they already represent the limits of the model's accuracy. 

The initial state for simulations will be X(l) = 1.0, X(2) = 1.0, in 
*For the constrained cases; ZC <_ 1.15, ZF 1.105, and T4 _< 1.08. 
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the less restricted quarter of the state space, so this control law will 
still be of interest. 

7.3 Model 3B Constrained 

As with model 3A, extremal controls appear throughout three quarters 
of the state space window. Again, X(l) = 1.0, X(2) =1.0 lies in the less 
restricted quarter. 

7.4 Model 4 Constrained 

Because different control space limits were used, this control law 
predicts better performance throughout the state space window. The ex- 
tremal controls do appear frequently, however, indicating that even better 
performance is possible with less restricted controls. Again, these control 
limits represent the range of accuracy of the model and may not be relaxed 
with confidence. 
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CHAPTER VIII 


CONTROLLER SIMULATIONS 


8.1 Introduction 

The performance of each control law generated was tested in two ways: 
by imposing the control law on the model from which it was generated, and by 
imposing it on the DYNGEN simulator. In each case, the initial state was 
the design point (u(l) = 1.0, u(2) = 1.0, x(l) = 1.0, x(2) = 1.0), and the 
final state was the control law target point (u(l) = 0.8, u(2) = 1.0, x(l) 

= 0.948434, x(2) = 0.928751). The simulation program used to impose the 
control laws on the models from which they were generated employs an Euler 
integration technique with a user controlled time step. This allows for 
the use of a smaller time increment In the steeper portions of a time re- 
sponse, without sacrificing the programming simplicity of the Euler tech- 
nique as compared to higher order numerical integration methods . 

The results of these simulations were not altogether satisfactory. 
Despite precautions taken to insure model accuracy, two of the model/control 
law systems are unstable. Experimentation with the integration time incre- 
ment, and the initial state, seems to indicate that the instability is not 
generated by the integration method. Application of the control laws to 
the DYNGEN simulator did not result in instability either, so this response 
is not inherent in the control law. 

This result does have a positive interpretation, however. Since it 
occured on the less complex models, 3A and 3B, one may conclude that the 
inclusion of more data points and derivative information, along with a more 
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sophisticated interpolation technique, does yield a better model. 


8.2 Model 3A with 3A-Controller 

The first application of a control law to its model resulted in insta- 
bility. After only 0.041 seconds, the states and outputs are beyond all 
physical limitations. Just prior to this overflow, u(l) oscillates sharply 
across a discontinuity in the control law, an effect Longenbaker [1] warns 
of in his work: 


"Interpolation will often lead to error when you are in- 
terpolating in a region of state space where the control 
laws change abruptly. Obviously, the optimal control 
which is desired is either one extreme or the other, and 
not something in-between." 

([1] 6.1 p. 78) 


This result would seem to indicate that model 3A, in its simplicity is in- 
capable of modelling the response of a system to a continuously varying con- 
trol. 


The states move outside the model’s range of accuracy, and the model 


breaks down almost immediately. 
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8 . 3 Model 3B with 3B-Controller 

This model/control law system is also unstable. After 0.027 seconds 
the states and outputs have exceeded physical limitations. In this instance, 
u(2) oscillates just before overflow, throwing the system into instability. 
Again, the indication is that the system is incapable of modelling the re- 
sponse to a continuously varying control , and model breakdown occurs al- 
most immediately. 
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8.4 Model 4 with 4- Controller 

This simulation yields a response that is stable. The states x(l) and 
x(2) converge to the target values after a small overshoot, and the outputs 
remain safely below constraint values. The states x(3) , x(4), and x(5) dis- 
play spikes and larger overshoots, but remain within acceptable limits. How- 
ever, the time estimate of 0.2642 seconds provided by the dynamic programming 
algorithm differs greatly from the actual target time of 0.5880 seconds. This 
can be attributed to the fact that the cost function is not linear, but its 
computation in the dynamic programming algorithm is by means of linear in- 
terpolation. Small errors in cost estimates near the target are both prop- 
agated and compounded outward through the state space. 
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8.5 Model 1 with 3A-Controller 

Despite the instability of model 3A with this controller, this simula- 
tion yields a smooth convergence to the target for states x(l) and x(2) . The 
controls display the same oscillatory behavior through a portion of the sim- 
ulation, hut the system remains stable. Analytically derived models embody 
physical relationships lacking in numerical models. This, and the complex 
coupling of a 16 state simulator provide for greater stability in DYNGEN 
than in model 3 A. 
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8.6 Model 1 with- 3B-Controller 

Again, DYNGEN is stable where the simpler model, 3B, is not. States 
x(l) and x(2) converge smoothly to the target, and output constraints are 
not violated. There seems to be little difference between the results from 
3A and 3B, indicating that the power law terra present in 3A is not particu- 
larly significant. 
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8.7 Model 1 with 4-Controller 

The response of the states and controls for this system are nearly iden- 
tical to those of model 4, the difference being that DYNGEN is about 0.05 
seconds faster. The outputs, however, vary greatly from model 4, indicating 
that the linear affine approximation of the outputs may be inadequate. De- 
spite this disagreement, the outputs do remain below the constraint values. 

As compared to models 3A and 3B, model 4 is certainly a step toward an auto- 
matically generated numerical model of a physical system. 
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CHAPTER IX 


SUMMARY AND CONCLUSIONS 

A dissatisfactory discovery of this work is that the accuracy of 
numerical models deteriorates so rapidly outside the range of their data 
points, that the state and control space limits must be confined to that 
range. The numerical model designer must be sure to include steady state 
and dynamic data from the entire area of state and control space that he 
wishes to explore. 

This rapid deterioration in accuracy is probably the cause of in- 
stability in models 3A and 3B. Having been derived from data at only 
two steady state points, these models are reliable only in the neighbor- 
hood of the line between these points . The transient generated by ap- 
plying the control law carries the state outside that neighborhood, and 
accuracy deteriorates to the point of total model breakdown. 

The accuracy and stability of the numerical models in this work were 
enhanced by the inclusion of more data, and a more sophisticated inter- 
polation scheme. However, this more complex scheme resulted in extrava- 
gant expenditure of cpu time. The cost of cpu time on Notre Dame's IBM 
370/168 is currently $275.00 per hour, while core allocation time costs 
$0.13 per K per hour. A typical run of 10 iterations of the dynamic pro- 
gramming algorithm with model 4 takes about 5 minutes of cpu time and 30 
minutes of real time. A reduction of cpu time of 20% to 4 minutes would 
result in a savings of $4.50. Assuming a similar reduction in real time 
to 24 minutes, this saving would- justify an increase in core usage to 88K 
k complete specification of 5 x 5 A matrices for 15 x 15 state space win- 
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dow, and 5 element g(u) function values for a 25 x 25 control space win- 
dow would occupy only 35K. This would result in a model in which there 
was no need whatsoever for interpolation, and probably a cpu time savings 
of greater than 20%. Values for eliminated states could easily be included, 
eliminating another major source of model inaccuracy. Pursuit of larger 
data bases, rather than more complex interpolation techniques seems more 
promising, in light of the availability and price of these two computer 
resources . 

The application of control laws to models involves interpolation be- 
tween adjacent control values. In the continuous portions of a control 
law, the difference interpolated is seldom more than a single control space 
increment. At a discontinuity in the control law, as stated by Longenbakei 
[1] and demonstrated in this work, interpolation is actually undesirable. 
Since control interpolation is either insignificant or counterproductive, 
it should be abandoned in future work in favor of a closest point scheme. 

A dynamic programming successive approximation algorithm was thor- 
oughly developed, programmed, and tested as part of this work. Designed 
to interface easily with any model, it should be of value in future re- 
search. 

In conclusion, the complex' interpolation scheme developed for this 
work did perform well, but the cost in programming and cpu time indicates 
that a shift to more data and less interpolation would be more profitable. 



APPENDIX A 


INPUTS FOR DYNGEN SIMULATOR 

This appendix includes the JCL used to run controller simulations on 
DYNGEN, the DISTRB subroutine that imposes the control law on DYNGEN, and 
the DYNGEN input data set, which includes the design point specification, 
three steady state requests, and a transient request. It is during the run 
of the transient that the control law is implemented. This JCL contains 
three separate jobs, one each for control law 3A, 3B, and 4. 
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//DYNGEN3A JOB ( CF ,F081 ). 7111530 15 , TIME = 1 0 * REG I0N=32UK » 
// N0TIFT=F907LB 
/* JOBRARM LXNES=20 


/❖ROUTE PRINT REM0TL2 
/❖ROUTE PUNCH REM0TE2 
/❖SETUP MU151 , RAT 

/❖SETUP PLOTS , NOCOUE 
//STLP1 EXEC PURTINIT 
//STEP2 EXEC PURTIN,TAPE=M015X 
//MOVEoSYSPRlNT DD DUMMY 
//STLP3 EXEC PURTSOUR 
//SOURCEoSY SPRINT UD DUMMY 

//SOURCEoSYSlN DU DSN=F9G7LB • DYNGEN, FORT s DISP=SHR 

//STEP4 EXEC PURTFORX t NAME=D JSTRB 

//FORT oSYSPRXNT DD DUMMY 

//LKED eSYSPRIN r DD DUMMY 

//STEPRA EXEC PGM=I EBGENER 

//SYSXN DD DUMMY 

//SY^PRXNT DD DUMMY 

//SYSUT2 DD OSN=X&CRDS»DISP=( NEW » PASS) , UNXT=DISK , SPACE= < 80 , (10,10) ) 
// DCB=(RECFM=FB,LRECL=60,BLKSIZE=80) 

//SYSUT1 DD USN=F9G7LB , DYNGEN »OATA ,OI$P = SHR 
//STEP5 EXEC PGM=IEBGENER 
//SYSIN DD DUMMY 
//SYSPRINT DD DUMMY 

//SYSUT1 DD QSN=&&CRDS,DlSP=(OLD,PASS) ,UNIT=DISK 
//SYSUT2 DD SYSOUT = A« f DCB=(RECFM=F»LRECL=80 «BLKSIZE=80 J . 

//STEPS EXEC PURTFGQ , NAME=DYNGEN1 

//G0«FT 0 7F0 01 N UD°DSN=2&PCH»UNIT = DlSK«SPACE=(TRKt (2*1) M DISP= ( * PASS ) 
//GO, FT08F0CIX UD D ISP= ( NEW , PASS ) 9 SPACE= ( TRK * ( 10 9 10 ) ) , UNIT=DISK 
//GOoFT09F001 DD * 

1/3A 
/ ❖ 

//GOoFTlOFOOl DD DSN=F9G7LB.DYN3AoDATA,DXSP=SHR 
//GO 9 FT11F0 0 i UD DSN-&&TRES,DXSP=( ,PASS> . 


J RINT DD DUMMY 

JT2 DD DSN=SaCROStDISP=(NEWtPASS)iUNIT=DlSK.SPACE=(80* 
DCB=(RECFM=FB ? LRECL=80fBLKSlZE=80) 

^ a » f~v > » 1 f-. n v/f.l/'r At n A Y nf OH-CUO 


(lotion 


//GOoFTlOFOOl DD DSN~F9G7LB»DYN3A , DATA,DXSP=SHR 
//GO 9 FT11F0 0 2. UD DSN=S&TRES , DX SP= ( « PASS ) , 

// UNXT=DISK*SPACE=(TRK, (10,2) ) 

//GOoSYSIN DD USN=&&CROS,DISP=< OLD, PASS) ,UNIT=DISK 
//TPLOT EXEC FORTHP 


//FQRToSYSPRINT 
//FORT. SYSIN DU 
//LKED, SYSPRINT 
//GO.SYSIN DD * 

0 0 0 9 0 0 8 

0 0 9 9 1 0 0 

0 o 9 « 1 e 0 

0 » 8 9 1 • 0 

0,85,1.1 

0 e 8 9 1 e 0 

0 . 7 « 1 » 0 

0 9 8 9 1 0 1 
0.85,1,1 
0o9,a,l 
0 0 85 , 1,1 
/❖ 

//GOoFTllFOOl DD 


DD DUMMY 

USN=F9G7LBoTFLOT,FORT,DISP=SHR 
OD DUMMY 


DSN=ZSTRES, DISP= ( OLD , DELETE ) 


//DYNGEN3B JOB ( CF, F081 ) ,711153015 , TIME=10 t REGI0N=32uK , 
// NOTIFY=F9G7LB 
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/❖JOBPARM LINES=20 
/❖ROUTE PRINT REM0TE2 
/❖ROUTE PUNCH REMOTE2 
/❖SETUP M0151,RAT 
/❖SETUP PLOTS, NOCOUE 
//STEP1 EXEC PURTINIT 
//STEP2 EXEC PURTIN,TAPE=MG151 
//MGVE.SYSPRINT OD DUMMY 
//STEP3 EXEC PURTSOUR 
//SOURCE. SYSPRXNT DO DUMMY 

/ /SOURCE » SYS I N DD USN=F9G7LB*,0YNGENoF0RT,DISP = SHR 

//STEP4 EXEC PURTFORX , NAME^DISTRB 

//FORT. SYSPRINT DD DUMMY 

//LKEO. SYSPRINT DD DUMMY 

//STEP4A EXEC PGM=lEBGENER 

//SYSIN DD DUMMY 

//SYSPRINT DD DUMMY 

//SYSUT2 DD DSN=&£CRDS,DISP=( NEW, PASS) , UNIT=DXSK, SPACE= { 80 , (10,10 > ) 
// DCB=(RECPM=FB,LRECL=80 ,BLKSIZE=80 ) 

//SYSUT1 DD USN=F9G7LB,0YNGEN.0ATA?DISP=SHR 
//STEP5 EXEC PGMslEBGENEK 
//SYSIN DD DUMMY 
//SYSPRINT DU DUMMY 

//SYSUT1 DD USN=&&CRDS, DISP= ( OLD » PASS ) , UNIT=DISK 
//SYSUT2 DD SYSOUT=A,DCB=(RECFM=F,LRECL=80,BLKSIZE=80) 

//STEP6 EXEC PURTFGO , NAME=DYNGEN1 
//LKED .SYSPRINT DD DUMMY 

//GO.FT07F001 DD DSN=&£PCH « UNI T=D ISK , SPACE= ( TRK , ( 2 , 1 ) ) ,DISP=< ,PASS) 
//GO.FT08F00I DD DISP= ( NEW « PASS ) , SPACE= ( TRK , (1 0 » 1 0 ) ) » UNI T=D I SK 
//GO.Fl 09F0G1 DD * 

1/3B 
/ ^ 

//GO.FTlOFOOl UD DSN=F9G7LB . DYN3B , DATA « OXSPrsSHR 
//GOcFTUFOOl DD DSN=££TRES, DISP= ( , PASS ] « 

// UNIT~DX$K«SPACE=(TRK, (10,2) ) 

//GQ.SYSIN DU DSN=&&CROS,DXSP=(OLD,PASS} ,UNIT=DXSK 
//TPLOT EXEC FORTHP 
//FORToSYSPRINT DD DUMMY 

//FORT. SYSIN DD DSN=F9G7LB . TPLOT „ FORT , DISP=SHR 
//LKED .SYSPRINT DD DUMMY 
//GO. SYSIN DD * 

0 , 0 , 0 ,6 
0 , 9 , 1 • 0 
0 e 9 , 1 e 0 

0 o 8 , 1 « 0 

0.85,1.1 
0 . 8 , 1 . 0 
0 , 7 , 1 o 0 
0 . 8 « 1 . 1 
0.85,1.1 
0 . 9 , 1 . 1 
0 o 85 , 1 o 1 
/* 

//G0.FTUFQQ1 DD DSN=SSTRES , DISP= ( OLD t DELETE ) 


//DYNGEN4 JOB ( CF * F081 >, 7111530 15 , TIME=10 , REGX0N=32UK , 
// N0TIFY=F907LB 
/❖JOBPARM LINES=20 
/❖ROUTE PRINT REM0TE2 


» 
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/❖ROUTE PUNCH REMOTE2 
/❖SETUP MO 151 » RAT 

/❖SETUP PLOTS, NOCOUE 
//STEP1 EXEC PURTINIT 
//STEP2 EXEC PURTIN, TAPE=M0151 
//MOVE # SYSPRXNT DD DUMMY 
//STEP3 EXEC PURTSOUR 
//SOURCE , SYSPRXNT UD DUMMY 

//SOURCE eSYSXN DD USN=F9G7LB#DYNGEN«F0RT,DISP=SHR 

//STEP4 EXEC PURTFQRX ,NAME=DISTRB 

//FORToSYSPRINT DD DUMMY 

//LKED sSYSPRlNT DO DUMMY 

//STEP4A EXEC PGM=IEBGENER 

//SYSIN DD DUMMY 

//SYSPRXNT DO DUMMY 

//SYSUT2 DD USN=&XCROS,DISP=( NEW, PASS) , UN XT=D I SK » SPACt = ( 80 , (10,10) ) , 
// DCB=(RLCPM=FB,LRECL=8D,BLKSIZE=80 ) 

//SYSUT1 DD QSN=F9G7LB,DYNGEN«DATA,DXSP=SHR 
//STEPS EXEC PGM=IEBGENER 
//SYSXN-DD DUMMY 
//SYSPRINT DD DUMMY 

//SYSUT1 DD USN-&&CRDS , DXSP= ( OLD , PASS ) , UNIT=DISK 
//SYSUT2 DD SYSOUT=A*DCB=(RECFM=F,LRECL=B0,BLKSI2E=80> 

//STEP 6 EXEC PURTFGO , NAME=0YNGEN1 

//GO^FT 0 7FOoi N UD U DSN=gSPCH,UNIT=DISK,SPACE=(TRK<! (2,1 ) ) «DISP=( ,PASS) 
//GO 0 FTO 8 FOOI DD D1SP=( NEW, PASS) ,SPACE=(TRK, (10,10) ) ,UNXT=DXSK 
/-/GO.FT09F001 OD * 

1/4 

/ 

//GO# FT10F001 DD DSN=F9G7LB 0 DYN4 0 DATA , QISP=SHR 
//GO# FT11F001 DD DSN=&STRES,DXSP= ( ,PASS ) , 

.// UNXT=OISK,SPACE=(TRK, (10,2) ) , 

//GOoSYSXN DD L)SN=&& CRUS , DISP= ( OLD « PASS > ,UNIT=DISI 
//TPLOT EXEC FORTHP 
//FORTOSYSPRINT DD DUMMY 

//FORT#SYSX'N DD DSN=F9G7LB ,TPLOT , FORT , DISP=SHR 
//LKED, SYSPRINT DD DUMMY 
//GO.SYSXN DD * 

0 0 0 « 0 08 

0,9, 1,0 
0,9, 1,0 
0 , 8 , 1 , 0 
0 # 8 5 , 1 « 1 

0 # 8 , 1 o 0 

0 0 7 , 1 # 0 
0 , 8 , 1 # 1 
0 , 85 » 1 # 1 

0 g 9 , 1 0 1 

0,85,1,1 

/ jjj 

//GO#FTUFOOl DD DSN = &STRES,DISP=(OLD, DELETE) 
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6 

O 


/ REPL NAME=DXSTRB 
/ NUMBER NEW1=1U.XNCR=1Q 
SUBROUTINE OISTRB 
IMPLICIT REALMS (A-H.O-Z) 

REAL*4 X ( 5 ) »U(2> .Y(3> ’ 

REAL+4 XMIN » XP5AX . XINC * XTAR 

REAl*4 vopt.uopt 

LOGICAL XNIT 

COMMON /COMALL/ COM (1062) 

COMMON XM IN ( 2 > »XMAX(2> • XXNC ( 2 ) « XTAR ( 2) 

COMMON VOPT( 15*15) rUOPT ( 15* 15*2) 

EQUIVALENCE ( TIME » COM C 993 ) ) « (DT*COM(994) ) . ( TF . CUM ( 995 ) ) 
EQUIVALENCE < WFB « COM (192 ) ) * ( A8 *,COM< 346 ) ) 

EQUIVALENCE ( XNHP * COM( 372 ) ) * ( XNLP . COM ( 374 ) ) 

EQUIVALENCE (P4«C0M(3B1> ) * < P7*COH( 389 ) ) 

EQUIVALENCE ( U4 , COM ( 382 ) > 

' EQUIVALENCE ( ZC , COM ( 30 0 ) ) « (2F« COM (136) ) . <T4?C0MU56) ) 
DATA TU/0*0/*X/5*l*0/tU/2#l,0/* Y/3*1.0/ 

DATA INI T/ « FALSE* / 

DATA XNHPN/0.118991D+05/ 

DATA XNLPN/G.987395D+04/ 

DATA P4N/ 0 » 2392990+02/ 

DATA P?N/0c235l42D+01/ 

DATA U4N/0*586468D+Q3/ 

DATA WFBN/2 » 75/ 

DATA A8N/2 9 9482558/ 

DATA ZCN/0.614303D+00/ 

DATA ZFN/0, 8333120+00/ 

DATA T4N/0»289208D+04/ 

IF(INIT) GO TO 10 
READ (9.202) MDL 
WRITE (11*202) MDL 
WRITEU1.200) TO. (X(I) *1=1*5) 

WRITE (11 . 201) (Ul I ) * 1=1*2) * (YU ) *1 = 1*3') 

READ ( 10 ) I7.VOPT.UOPT 

READ ( 10 ) XMINi XMAX* XINC » XTAR 

INIT“ • IRUE « 

10 CONTINUE 

X(X>=XNHP/XNHPN 

X(2>=XNLP/XNLPN 

X ( 3 ) =P4/P4N 

X(4}=P7/P7N 

X(5)=U4/U4N 

CALL CNTRL(X.UfV) 

WFB=U(1)*WFBN 

A8=U(2)*A8N 

Yll)=ZC/ZCN 

Y(2)=ZF/2FN 

Y(3)=T4/14N 

WRITE UX *200) TINE. (X(I) *1=1*5) 

WRITE (11*201) <U(X) *1=1*2) *(Y(I)* 1=1*3) 

RETURN 

200 FQRMAT(F6„5.b£14 0 6> 

201 FORMAT ( 6X . E5E14 »6 ) 

202 FORMAT (A4) 

END 

SUBROUTINE CNTRL(X.U.V) 


C 

C 


THIS ROUTINE INTERPOLATES CONTROLS FROM THE 



ENTRIES OF THE OPTIMAL CONTROL LAW. 
REAL X<5) tU(2) 

COMMON XMIN (2) ? XMAX(2) »XXNC<2) sXTAK(2) 
COMMON V0PT(15»15) »UOPT U5.15«2) 

NX = 15 

FR1H= (X(1)»XMIN{1) )/XINC(l)^loO 
IX1L=IPIX(FR1H) 

FRXH=FR1H-FL0AT< XX1L) 

IF UXlL.GT.NX) IX1L=NX 
IXXH=IX1L+1 

IF ( IXlH.GT .NX) IX1H=NX 
FR1L = 1 » 0 “FR1H 

FR2H=tX(2)-XMIN(2) ) /XINC ( 2 ) • 0 
IX2L=IFIX(FR2H) 

FR2H=FK2H-FL0AT(IX2L) 

IFUX2L.GT.NX) IX2L=NX 
IX2H-IX2L+1 

IFUX2H.GT.NX) XX2H=NX 
FR2L = 1 • 0-FR2H 

V=FRlL*FK2L*V0PT ( IX1LUX2L) 
Utl)=FRlL*FR2L*UOPT ( IXXL * IX2L ».l ) 
U(2)=FrtlL*FR2L*UOPT(IXlL*IX2Lt2) 
VcV+FRlH*FK2L*V0PT( IX1H* IX2L) 
U<l)=U<l>+FRlH*FR2L*U0PTUXlHtIX2Ltl) 
U<2)=U(2>+FRlH*FK2L*U0PT(lXlHt IX2L*2) 
V=V+FRXL*FR2H*VQPT ( IX1LUX2H) 

U{ 1 ) su ( 1 ) ■s-FR1L ! < c FR 2H*UQPT (IXlL*IX2Htl) 

U<2)=U(2)+FR1L*FK2H*U0PT(IX1L»IX2H«2) 

V=V+FR1H*FR2H*V0HT( IXiHt IX2H) 

U(1>=U(1)+FR1H*FR2H#U0PT(IX1H«IX2H*1) 

U(2)=U(2)+FR1H*FR2H*U0PT( IXIHt XX2H»2) 

RETURN 

END 
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WFB 

A8 

XNHP 

XNLP 

PA 

P7 

U4 

ZC 

ZF 

TA 

FG 

THEEND 

SDATAIN ISP00L=2»FAN=*TRU£, «SI=, FALSE, » IDES=1 » HODE=0 • IDUHP=1 , IAHTP=0 
IGASMX=2 t ITKY8=20U,FXFN2M=, FALSE* »FXM2CP=. FALSE , *AFTFAN=, FALSE • t 
DUMSPL=. TRUE, *T0LALL=l,E-6 tDELFG=X* 0 *DELFNs:l,0 «DELSFt = 1,0 t 
PCNF0S=1U2,31* 

PRFDS=2. 99b «ETAFDS=. 8499 »PCNCDS=98, 73 tPRCDS=B. 462 iETaCDSs, 8136 • 

T40S = 2892c OA, 

ETABDS=l.QQ«DPCODS=, 0561 »ETHPDS=, 8713 «ETLPDS=, 9021 *DPDUOS=. 0584 t 
T7DS~3583, 6 « 

ETAADS= , 84-30 , UpAFUS= , 0599 1 AM55= , 283 « AM6= , 243 1 CVFiNOZ= • 9494 i 
WAFCDS=221,373, 

.UACCOS=54,988«HPEXT=0,0« AH=0,0« AL.TPsO.Ot PCBLF=0,0«PCULC=*18* 
PC8L0U=,208t 

PCBL0B=0, 0 » PCBLHP= , 726, PCBLLP=, 066 *AM23=.l 70 *ZFDS=o 8*33, ZCDS=. 8143* 
TFHPDS=50 » 0 * CNHPDS=2 , 0 » TFLPDS=130 , 0 » CNLPDS=2, 3 « XNHPDS=10070 , » 
XNLPDS=9651* » 

P«IHP=3*80tPMILP»4,50*VFAN=2,31,VC0HP=1.65«VCOf1B=l*6o«VHPTRB=,505* 
VLPTRB=,blU» VAFTBNS49*77» VFUUCT=10,08. SEND 

SDATAIN IDES=Q,INIT=1* ITRAN=0 * M0DE=2 , 

A8=2 9 9482B58.tWF6=2.75 SEND 

SDATAIN IDES=O»lNlT=l*ITRAN=0.nODEs2» 

A6S2.1442* MFB=2,0 SEND 

SDATAIN IDESsO*INIT=:lt ITRAN=O.HODEs2» 

A8=2,9187732iHFB=2,3375 SEND 

SDATAIN M00E=2»ITRAN=:1«INIT=1* 

OT=0.01«QTPKNI=0.1*TF=0,8i 
PMIHP= 3e80sPMILP~A»50 i 
VAFTBN=8,0,VCUHB=l,0.VCOHPal.0iVFAN=l,0t 
VFDUCT=2.0« VHPTKB=l*Oi VLPTRBsl.Oi 
XNHPDS=10070»U,XNLPUS=9651,0« 

WFB=2»75,A8=2*9A82558 SEND 

/* 



APPENDIX B 


MODEL 3 AX GENERATOR PROGRAM 

This appendix includes the Model 3 program, and the two input data sets 
used to generate derivative values for Models 3A and 3B. This program is of 
a form that allows it to be compiled with either the dynamic programming al- 
gorithm, or the general system simulator. 
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SUBROUTINE INIT(MDL) 


100 

200 

201 

202 

203 

204 

205 

206 


INITIALIZE ALL VALUES NEEDED 
BY THIS MODEL* 

REAL ALPHA,BETA,XD(5) * XW ( 5 ) * AD O * 5 ) * AW< 5 * 5 ) , C ( 5 , 6 ) 

COMMON/ MDL/ ALPHA, BETA , XD , XW , AD , AW , C 

READ (1,100) MDL 

WHITE ( 6 * 200 ) MOL 

READ (1*3) ALPHAtBETA 

WRITE ( 6 ,201 ) ALPHAtBETA 

READdt*) XD 

WRITE (6,202) XQ 

READ (It*) XW 

WRI TE ( 6 , 203 ) XW 

READ (It*) ( (AD(ItJ) ,J=1«5) , I=l»b) 

WRITE (6, 204) ( ( AD( I ,U) t J=1 ,5) t 1=1,5) 

READdt*) ( ( AW ( 1 1 J ) * J=1 15) , 1 = 1 » b ) 

WRITE (6 ,205) ( ( AW( I * J ) , J-i , 5) ,1 = 1,5) 

READ (It*) ((C(XtJ)t vi=l 1 6 ) * 1=1 » 5 ) 

WRITE (fa,206) ( (C(ItU) t J=lt6) ,1=1,5) 

RETURN 
FORMAT ( 1A4) 

FORMAT { IX , * DATA FOR MODEL *,A4) 

FORMAT (IX, * ALPHAr* ,E16.6, * BET A= * ,E16*6) 

FORMAT ( IX, 'XD=* ,T20 ,5E16,6) 

FORMAT ( IX, «XW=* ,T2U*5E16 0 6) 

FORMAT (/, IX, 1 AD= 5 ,5(T20,5El6o6,X) ) 

FORMAT ( /,1X, ’AW=» ,5 (T20,5E16,6d ) ) 

FORMAT (/, IX, f C=' ,5(T20,6E16*6t/) , «1» ) 

END 


SUBROUTINE STDST(U,X) 

THIS ROUTINE EVALUATES G(U), 

THE STLADY STATE FUNCTION * 

REAL U ( 2 ) # X ( 5 ) , US ( 2 ) 

REAL ALPHA,BETA,XD(5) ,XW(5) ,AD(b,5),AW(5,5),C(5,6) 
COMMON/MDL/ ALPHA , BETA , XD , XW , AO , AW , C 
DO 10 1=1,2 
US ( I ) =ALPHA*U ( I ) +BETA 
10 CONTINUE 
DO 20 1=1,5 

X< r>=C(X»D*US(l)+C(I,2)*US(2) 
8-e>C(I,5)*US(l)**C(I,3)*US(2)**C( 1,4) 

( 1 , 6 ) 

20 CONTINUE 
RETURN 
END 
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SUBROUTINE XDOT(X,U»DX) ' 


THIS ROUTINE EVALUATES 
THE DERIVATIVE FOR ANY 
STATE AND CONTROLS* 


10 


REAL X(5)iU(2>'0X(5)*G(5)«UK(5)«USy(2)«A(5«5) 

REAL ALPHA t BETA , XD ( 5 ) , XW ( 5 > iADOf5)*AU(S«5)tC(5«6) 
COWMON/MDL/ ALPHA , BETA , XD , XW » AD . AW , C 
CALL STDST ( U « G ) 

DO 10 J=l*5 
WK{ J)=X{ J>-G( J) 

CONTINUE 

CALL APIAT(X»U« A) _ „ _ TrM 

CALL VMULFF(A,WK»5t5,lt5t5,DX«DoXER) 

RETURN 

END 


SUBROUTINE AMAT ( X » U « A ) 


THIS ROUTINE GENERATES THE 
A MATRIX FOR ANY STATE 
AND CONTROLS* 


10 


REAL X(5),U(2),A(5«5) 

REAL ALPHA, BETA ,XD<5> , XU ( 5 ) t AD < 5 , 5 ) ♦ AW ( 5 , 5 ) , C ( 5 » 6 ) 
COMMON/MDL/ ALPHA , BETA , XD , XW » AD , AW , C 
DO 10 J=1 , 5 


FRW=(XO(J)-X(U) )/(XD( J)-XW( J) ) 
FRD= ( X ( J ) «*XW ( J ) )/(XD( J)-XU( J) ) 
DO 10 1=1,5 


A ( I » J ) =AW ( I * J ) #FRW+AD ( I , J ) #FRD 

CONTINUE 

RETURN 

END 
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3A 

1,1, -Ool 

1 . 0 . 1 , u 1 1 . 0 , 1 , U , 1 , 0 

0 ,9,0.789710.7381,0.9401*0.9454 

-3 .8,-1.277,2,067,-1.152,1.448 

2.748, -5,39. 1.585, -1.991 i 1.071 

377 , 9, 49. bl, -264.9,66. 807, 78, 91 

31. 2b, 139. 39 *-6.269, - 88. 69-, 27. 83 

-176. 5, 23. 91, -10.27,-37.4,-246.7 

-4.744, -1.3888 ,3. 2468, -1.4591, 1,1969 

0. 82186, -2b, 726, 2 .5585, -1,8609, 0,45548 

475.73,137.55,-328,91,27,791,91,495 

-50. 103, 110 .91 ,63 .186, -116 .69, 6 .2883 

-186,77,-67.682,-41,681,24.586,-243.23 

0,24267, -0.0 0218, 1.90082, 8. 09916, 0.02864, 0.73088 

1.01593.0. 85407.0.89872.0.66919,-0,81679,-0.05121 
0.73445,0,10133,6.90586,3,09409,0,011495,0,15272 
0.77234,-0,35905,2,49867,2,87415,-0,075198,0,66191 
0,39503 ,-0,27262,-3, 44682, 13, 4468, 0,01838,0,85921 


3B 

2.31778,-1.31778 

0 , 9 , 0 , 7897 , 0 , 7581 , . 0 , 9401 1 0.9454 

i9 « 1 • 5 85 » -1 . 991 , 1 . 071 

»17l 6 A 1 ?l , u? f “^ 26 l’“ Q8e69,27 « 8 3 
'*176,5, 23, 91, -10, 27, - 37. 4, - ^46. 7 

0 4 8?^rtfi" ;l p2 8 7fA 3 o 2 -Ao , '' :L * 45yl,1 ' 1969 

4§ 2 J ^ 2b * 728 * 2 • 0585, -1 , 8609,0.45548 

-40 e in4* L ?Tn‘ 3 qi'*'f 2b i^a' 27<,79l,9 ^ et|9 ' 3 

n?Z?£T 7 !“«T* ^ 82 ’”41,681,24.586,-243.23 

^’5351 ,.-0.1208, 1.0,1. 0*0, 0,0. 5857 
n * “2 * 1,0,1,0,0,0,0.9053 

0,2962 ,-0.2099, 1,0,1,0,0.0,0,9157 



APPENDIX C 


MODEL 4 AX GENERATOR PROGRAM 

This appendix includes the Model 4 program, and the input data set used 
to generate derivative values for Model 4. This program is of a form that 
allows it to be compiled with either the dynamic programming algorithm, or 
the general system simulator. 
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SUBROUTINE INIT(HDL) 

INITIALIZE ALL VALUES NEEDED 
BY THIS MODEL, 

INTEGER ORD ,N 1 1 A ( 5 ) 

REAL U2 , U1 ( 5 > , XI ( 5 « 5 ) « DXI ( 5 , 5 ) 

REAL X(5) »U(2> tAl(5»-5)«A3(5*b) 

REAL AP (5,5) » AM ( 5 i 5 ) 

COMMON/MOL/ 0RD,N,XA«U2,UI,XI, 

COMHON/MDL/ A1.A3.A5 

COMMON/MOL/ AP , AM 

DATA MDL^/ * 4 V 

MUL=MDL4 

READU,*) ORD.N 

READ ( 1 ? * ) U2 

DO 10 1=1, N 

READ( Is*) UKI) 

■READU**) (XKl, J) *J=l«ORO) 
READ {It*) (DXHI.J) ,U=l»ORD> 
REAO(l«*> I A { I } » DEL { l ) 

IFUA(X) , NE , 1 , AND » I A ( I ) ,NE,2) 
READ ( 1 . * ) (XIP(I»J),J=l,ORD> 

IF ( XA ( X ) «EQol ) GO 10 10 
READ ( It*) (XIM(IiU) , U=1.0RD> 

10 CONTINUE 

READ ( 1 1 100 ) ( ( Alt I, J) , J=l,ORD> 
READ(lilOO) ( (A3U.J) ,J=l,ORD) 
READ 11,100) { (A5( I, J) ,J=l,ORD) 
READ (1*100) ( lAPlItd) » J=l,ORO) 
READ (1*100) ( (AM(I.J) »U=X»ORD) 
RETURN 

100 FORMAT ( 5£X3,6} 

END 


THE NEXT EIGHT ROUTINES ARE THE 
HERMITE INTERPOLATION METHOD OF 
HILDEBRAND, MODIFIED AS PER THE TEXT. 

REAL FUNCTION L(I,U) 

INTEGER ORD « N , I A ( 5 ) 

REAL U2,Ul(5),XI(5,5),DXI(5,5),xlP(5,5),XIH(5,5) .DEL (5) 
REAL X ( 5 ) , U ( 2 ) « A1 ( 5 , 5 ) , A3 ( 5. 5 ) , A5 ( 5 , 5 ) , A ( b « 5 ) 

REAL AP(5,5)tAM(5,5) 

COMMON/ MDL/ ORD » N , I A » U£ , U1 , XI , D *1 , XIP , XIM . DEL 
COMMON/MDL/ A1.A3.A5 
COMMON/MDL/ AP , AM 
L=1 o 0 

DO 10 K=1»N 

IF { K , EQ , I ) GO TO 10 

L=L*(U(1)-U1(K) )/(Ul(I)-UKK) ) 

10 CONTINUE 
RETURN 
END 


,XIP(5«5) ,XIM(5»5) , DEL ( 5 ) 
»Ab(5,5) ,A{5,5> 

DXl , XIP , XIM, DEL 


SlOP 16 


, 1=1. ORD) 
» 1=1, ORD) 
. 1=1, ORD) 
, 1=1 , ORD) 
, 1=1, ORD) 



20 

10 


REAL FUNCTION DL(X«UU> 

KEAL*^U2 *U1?M « XI [lt5) »DXI ( 5« 5) » XIPt 5 »5> I f *DEL (5) 

REAL X ( 5 ) 4 U ( 2 ) ?Al(5»b ) 4 A3 l 5 4 b ) 4 Ab ( 5 4 5) « A { b «5 ) 

COMMON /mBl/ OROsNslAt U2 tUl*Xi»DAl?XIPvXlM» DEL 
COMMON/MDL/ A1»A6«A5 
COMMON/MDL/ AP « AM 
DL=0 » 0 

00 10 K= 1 *N _ , (1 

IF(K.EQ.I) GO TO 10 

PRD=l e 0 

00 20 KK=1,N ^ 

IF ( KK « EQ * I ) GO TO 20 
TERM=UU-U1(KK) 

IF ( KK « EQ « K ) TERM=1*0 

PRD = PRU^TERM/(U1.(I)-U1(KK) ) 

CONTINUE 

□L=OL-i“PRD 

CONTINUE 

RETURN 

END 
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REAL FUNCTION R(I*U) 

RE AL^U2 eUl?5^«XI(l!5) »DXI ( 5» 5 ) » XlP ( 5*5 ) *XIM ( 5 *S) *DEL(5) 
REAL X(5) *U(2) 4Al(595)«A3(5«5)*Ab(5*5) , iA(5s5) 

COMMON/MDL/* 6rd!n»IA«U2*U1iXX*DxI ,XIP*XIM*DEL 
COMMON /MOL/ Al4A3?A5 

COMMON/MOL/ AP.AM 

'R=l.O-2*0#DL( J*U1(I) )*(U(l)»UUi) ) 

RETURN 

END 


REAL FUNCTION SU 4 U) 

INTEGER ORD » N « I A ( 5 ) VT „,, K nc1 

REAL U2 ? U1 { 5 ) « X I ( b i 5 ) * DXI < 5 » 5 ) » X I P ( 5 ♦ 5 ) *XIM(5»b) »DEL<5> 
•REAL X(5)*U(2)*All5*5)tA3(5*b>»Ab(5i5)«A(5«5) 

COMMON/MOL/ ORD t N » I A » U2 » U1 « XI *0*1 « XIPi XIM* DEL 
COMMON/MOL/ A1«A3,A5 
COMMON/MDL/ AP« AM 
S=(U(1>-U1(I> ) 

return - 

END 


REAL FUNCTION HCI»U) 

REAL G U2,Ul?5^txi!ll5) , DXI < 5 * 5 > . XlP { 5 *5.) 4XIM<5»5> ,DEL*<5) 
REAL X(5)«U(2>iAl(5f5}<A3(5*b)tAb{5*5)»A(555) 

REAL AP ( 5 » 5 ) * AM (5 1 5 ) 

REAL- L 

COMMON/MDL/ ORD * N 4 1 A 4 U2 1 U1 * XI »DxI « XIP t XIM * DEL 

COMMON/MDL/ Al4A3*A5 

COMMON/MDL/ AP 1 AM 

H=R(ItU)*(L(IiU))**2 

RETURN 

END 
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REAL FUNCTION HBUtU) 

INTEGER ORD ? N 1 1 A ( 5 ) = c , 

REAL U2 » U1 ( 5 ) « X I ( b * 5 } t DXI l 5 ? 5 } , xip ( 5 , 5 ) t X IM { 5 « 5 ) tDEL(5) 
REAL X(5)tU(2)*Al(bt5)«A3(5*5)»Ab(5e5)tA(5*5) 

REAL AP(5»5) »AM(5» 5) 

REAL L 

COMMON/MOL/ ORD »N » I A » U2, U1 « X I * Dx I « X I Pc XIM * DEL 

COMMON/MDL/ A1,A3»A5 

COMMON/MDL/ AP,AM 

HB=S(ItU)*(L(I*U) )**2 

RETURN 

END 


10 


REAL FUNCTION BOX(ItJtU) 

INTEGER 0RD»NiIA(5) 

REAL U2 » U1 £ 5 ) « X I ( 5 ♦ 5 ) . DXI < 5 « 5 ) » xip £ 5 1 5 ) t XIM £ 5 t 5 ) ,DEU5) 

REAL X(5)*U(2)«Al(5«5)fA3(5tb)«A5(5'5)tA(5t5) 

REAL AP ( 5 » 5 ) » AM ( 5 » b ) 

COMMON/MDL/ QRDtNcIAcU2tUl,XItDxI*XIPtXIM» DEL 

COMMON/MDL/ Alt A3 » A5 

COMMON/MDL/ AP t AM 

IF{ I A £ I ) ,£Q , 2 ) GO TO 10 

22£r^, 0+{U<2, " U2) * (XIP(I » U)-XIUt J) )/(DEUI)*XI< It J) ) 

KE.T UKN 
CONTINUE 

o?R?^ < ’ 1 0 1 S{ U( £ ) 7 U2, * (XIP{I ' J >~ ximi ' J}) ''< 2 <>0*DE*-(I)*Xl<ItJ)K 

&<U<2>-U2>**2*<XIP<I»J>+XIM<ItJ)-2.0*XI(ItJ>)/ 

& ( 2, Q#DEL ( I ) #*2*XI ( I « J ) ) 

RETURN 

END 


REAL FUNCTION Y(JtU) 

INTEGER ORD « N 1 1 A { 5 ) 

REAL U2.U-K5) »XI<5t5) tDXI<5»5) *XIP(5 i 5) *XIH(5«5) cDEL£5> 
REAL X(5) ?U(2) « A1 ( 5 « 5 ) t A3 £ 5 » 5 ) t Ab ( 5 » 5 ) c A ( 5 c 5 ) 

REAL AP(5*5) » AM£5.5) 

COMMON/MDL/ ORD , N * I A i U2 t U1 t XI 9 Dxl t XIP « XIM t DEL 
COMMON/MDL/ A1»A3»A5 
COMMON/MDL/ AP,AM 
Y=0c 0 

DO 10 K=ltN 

Y = Y-s-H£K,U)*XI (K» J)*BOX(Kc J*U)+HB<KtU)*DXI£Kf J) 

10 CONTINUE 
RETURN 
END 


10 


SUBROUTINE STDST(UtX) 

THIS ROUTINE EVALUATES G(U). 

THE STEADY STATE FUNCTION, 

INTEGER ORDtNtlA(S) 

REAL «DEL(5) 

coBSoK/mol/ K?;3:i8 ,U2,ul,xl,B,a,XIP,xl " ,OEt 

COMMON/MDL/ AP , AM 
DO 10 J=l ,5 
X( J)=Y( JtU) 

CONTINUE 

RETURN 

END 
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5 

10 


SUBROUTINE XDQT(X,U,DX) 

THIS ROUTINE EVALUATES THE 
DERIVATIVE FOR ANY 
STATE ANQ CONTROLS. 

INTEGER ORD ,N,IA(5) „ rr ,, c VTM , C 

REAL U2*U1(5) ,XI<5,5) tDXI(5«b) tXiP(5*5) »XIM(5 
REAL Xl5>,U(2),Al(S,b),A3<5,b) } Ab{5,b>«A{b,b) 
REAL AP{5,5) * AM ( 5 » 5 ) 

REAL DX ( 5 } « G ( 5 1 , WK ( 5 ) , USV ( 2 ) VTO vTM 

COMMON/ MU L/ OKD , N , 1 A t U2 , U1 , XI t DX I » XIP * XIM ? DEL 
COMMON/MUL/ A1*A3»A5 
COMMON/MDL/ AP « AM 

DATA USV/2*1.Q/,G/5*1,0/ , , t 

XF ( U ( 1 ) oEQoUSV(l) .AND.U(2) .EQ.USV(2) ) GO TO 5 


USVU)=U(1) 

USV ( 2 ) =U ( 2 ) 

CALL STDST(UiG) 

CONTINUE 
DO 10 0 = 1,5 
WK{ J>=X(U)~G{ J) 

CONTINUE 

CALL AMAT ( X , U * A ) 

CALL VMULFF (A,WK»5*5il«5,5«DX,5,IER) 


RETURN 

END 


»5)i DEL ( 5 ) 


SUBROUTINE AMAT(X.U.A) 

THIS ROUTINE GENERATES THE 
A MATRIX FOR ANY STATE 
AND CONTROLS. 

INTEGER ORD* N, XA (5 ) 

REAL U2»U1(5> *XI(5»5> »DXI<5«5> fXlP(5*5) * XIM ( 5 « 5 ) * DEL { 5 > 
REAL X(5).U(2)»Al{5»5}.A3(5»5).ft5{5«5)*A{5*5) 

REAL AP (5.5) »AM(5,b> 

COMMON/MDL/ ORD * N i X A , U2, U1 , X I , Oxl , XIP * XIM » OtL 
COMMON/MDL/ A1»A3.A5 
COMMON/MDL/ AP. AM 

FR1=(X(1)-XX<3, 1) ) * ( X ( 1 ) ~XI ( 5 , 1 J ) 
&/({XI(ltl>-XI(3,l))*(XI(l,l)-XXl5tl))) 

FR3s(X{l)-XHlil) >*(X(l)-XI(b,ll) 
a/{ (xi(3ti)-xi(i«i) )*(xi(3»i)-xiis»i>n 
FR5=IX(-1>-Xni*li >*(X(1>-XI(3.1M 
&/( (XI (5,1) -XI (1,1) J*(XH5«l)-XIt3.1) ) ) 

UlP=U2-i-DEL( 1 ) 

U1M=U2°0EL(1) 

FRP=(U(2>-U2)*(U(2)»UIM)/( (U1P-U2) * (U1P-U1M) > 
FR=(Ut2)*>UlP)*(U<2)-UlM)/{ ( U2-U1P >* ( U2-U1M >) 
FRM=(U{2)-'UlP>*(U(2)-U2>/( ( UlM-UiP ) * ( U1M-U2 ) ) 

DO 10 1=1, OKD 
DO 10 U=1 , ORD 

AU,U)=A1(I, J)*FR1+A3( I,J)*FR3+a5(I, J)*FR5 

A(I»J)=A<I,J)*< FRP*AP ( I , J ) /A1 < I , J > +FR+FRM*AM ( I »J)/Al( I tU) ) 
10 CONTINUE 
RETURN 
END 
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5*5 

1,0 

1,0 

0,10000 OE+0 1 0,100000E+01 0,100000E+Ql OelOOOOOE-f-Ol UoXOOOOOE+Ql 
0 , 258696E-J-00 0.232QQ6E+0Q 0,741347£+00 Q,536700E+Q0 u,368200E+00 

2,0, 13567 

0 , 102201£^01 0 , 1 03699E-5- 0 1 Q,1Q0909E+01 0 ,922052E+00 U,X001X5E+01 
0,983596E+00 0 , 944368E+Q0 0,972825£+0Q 0,106272£+01 u,102339E+01 
0,945455 

0 « 984789E+00 0 o 985735E*00 0 »9566Q4E-(-00 0 , 96B504E+00 U.98X774E+00 
0,280586E+00 0,266837E+00 0.798173E+00 Q.582468E+00 u , 338984E+00 

2,0,13567 

0 , 1010Q9E+01 0 « 102221E+01 0 » 969378E+0 0 0 , 894319E+G 0 U,980202E + 00 
0,968989E+00 Q,927607E + 0G 0.931285E+G0 0 , 1 02997E+01 U,100344£+01 
0,872727 

0 , 966384E+00 0 , 959990E-5-00 0 , 898753E+0 0 0,923883E + 00 U,956888E+00 
0 ,23490 1E+OU 0,417786E+00 0,794424E+00 0 „ 635940E+0Q U „ 345763E + 0 0 

2,0,13567 

0 o 990008E+00 U „ X0002QE+01 0.913027E + 00 0 , 855471E+00 U.953457E+00 
0,948022E+00 0,9039l2E + 00 0 « 874094E-f-Q0 0 , 984805E+00 U,976756E + 00 
0,8 

0 , 948434E+0 0 0 ,928751E+00 0 ,840108E+00 0,876949E+00 U.931145E + 00 
0 ,265795E+0O 0 , 442441E+0 0 0 ,827564£+00 U , 656344E+00 u .353039E-S-00 

2,0,13567 

0 « 967771E-5-Q 0 0 , 975912E+00 0 , 853931E+Q0 0»814828E+0Q U.926353E+00 
0 , 926608E+0 0 0 , 676112E + 0 0 0 , 8162XGE+Q0 0 , 938098E-8-0 0 u,948750E + 00 
0,727273 

0,928298E+00 0 , 895714E+G 0 0 ,779278E+00 0 , 828553E+00 U , 904594E+0 0 
0«286381E+00 0 »464X84E*00 0,843807E+00 0 , 673D48E+ 0 0 0,375500E + 00 

2,0,13567 

0,94792lE+00 0 , 940920E+Q0 0,794242E+00 0,770606E+00 U,898666E+00 
0 e 907843E+0 0 0 , 845919E+00 0 ,758233E + 00 0 , 887 0 10E*00 u,919949E+00 
”0o439400E+01“0, 988129 E+ 00 0 ,252252E<-01~0 , 157586E+01 U,123698E+01 
0 E 940742E+00«»0 , 60930 0E+01 0,266571E+0X~0 , 184 044E+01 U ,707764E + 00 
0 « 506886E*s-05 0 , 519848E-i-02~ 0 , 350 590E + 03 0 , 146869E+ 0 3 U,116621E+Q3 
“0 o 585698E-S- Q2 0 . 128300E + 03 0 .607624E+02-Q . 121120E+03 U . 113162E+02 
®0 « 159711E+05 0 ,26b56bE + 02*'0 »189252 Ex 02»'0 ,446012E+02“U ,242670E + u3 
“0,471740E'S-0l”0,124636E-}-01 0 • 298608E+01-0 , 187860E+01 U . 143928E+ 0 1 
0 ,5526 0 1E + 0Q~0 ,35739 0E+ 01 0 ,2782 04E+Q1-Q ,1721 09E + 01 u „ 723029E-S- 0 0 
0o450922E + 03 0 , 813586E+02-0 , 332560E+03 0 , 130U44E-S-03 U,113546E + 03 
®0 o 517156E + 02 Q,S88692E + 02 0 , 585464E+02-0 , 119240E+03 U , 102770E + 02 
-0,130919E+0 3-0 o 451477E+01-0,341646E+02-0,372193E + 02*-U,247470E + 03 
•»0p478520E-{-01“0,121176E + 01 0 ,341043E+01-0 ,209109E+0X u,135367E + 01 
0,499094E+00*»0, 347050E+01 0 , 258839E-S-0 1 = 0 , 140777E+01 U, 625560 E + QO 
0 o 398468E+ 03 0 , 805642E+02- 0 , 318460E+03 0,113763E+03 U,1021X3E+03 
«0 c 434077E+02 0,8442265+02 0 , 561960E+02-0 • 116470E+03 U , 944216E+01 
•=OnX13790E+03°0, 53623 QE+ 01-0, 455125 E+G2-U, 323419E + 02 ~u «242430£+03 
» 0, 73132 OE+0 1=0, 635114E+00 0 , 284068E+ 0 1 -0 , 166920E+ 01 u ,150321£+01 
Q,X5Q544E+01»0,549860E+G1 0 ,28101 5E + G 1 = 0 , 227371E+01 U,753975E+00 
0,526328E + 03 0 , 521252E+02-0 « 329620E+Q3 Q,121962E+03 U,127924E + Q3 
«0,61206l£+Q2 0 , 109032E+03 0 ,571358£+02=Q ,X223-70£ + 03 U,991496E + 01 
**G,X20121E+O3 0, 167193E+02~G, 340 961E+ 02=0 , 332664E+02-U , 247680 E+03 
-0,324940E+01»0,914027E+00 0 « 194351E+01-O ,139379£+01 U.147637E+01 
0„10Q929E+QX»Q „360220E+01 0 , 259057E+01~0 , 227073E+0 1 U,665059E400 
0 « 566953E+Q5 0 , 152253E+03-Q « 54Q260E+Q5 0,3l5925E + 02 u,116803E + 05 
»0« 59 54 59 E -8- 02 0,X19995E+03 0 ,63552X£*Q2~0 , 123660E + 03 u,102124E+02 
-0, 22926 9E+03«0, 61997 3E + 02 ra O,938650E+01 0 , 19816 0E + O2-U ,253130E-t-03 



APPENDIX D 


GENERAL SYSTEM SIMULATOR 

This appendix includes the JCL used to run controller simulations on 
the model from which the control law was generated, and the general system 
simulator program, TRES. This program uses a Euler integration technique 
with a user specified time increment. A data set of time response values 
is generated, with each record consisting of time, state, control, and 'out- 
put values. This JCL contains three jobs, one each for models 3A, 3B and 4. 
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//TPLQT3A JOB ( CF , F081 ) t 711133015 . N0TIFY=F9G7LB , 

// time=i 

/*rqute print remoter 
/&SETUP PLOTS i NOCODE 

//TRES EXEC FORTH, LI84=* IHSL. SINGLE* ■ 

//FORT, SYSPRINT OD DUMMY 

/ /FORT .SYSIN DU DSN=F9G7L8 . TRES, FORT , DISP=SHR 
// OU OSN=F9G7LB 0 MOOEL3 # FORT(DiSP=SHR 

//LKED, SYSPRINT DO DUMMY 

//GOoFTOlFQQl UD DSN-F9G7L8 , M0DEL3A , DATA , DISP=SHR 
//GO. SYSIN OU * 

1.0. 1.0.0.005 

41.0. 001 
/* 


//GOeFTlOFOOl DD DSN=F9G7LB,DYN3A.DATA«DISP=SHR 
//GO.FTllFOOl 00 DSN=SXTRES,DISP=( .PASS) » 

// UNIT=OISK»SPACE=(TRK«UO»2) ) 

//LIST EXEC T SOLI ST , NAME= 9 &&7RES * , COND=EVEN 
//TPLOT EXEC FORTHP 
//FORT, SYSPRINT 00 DUMMY 

//FORT. SYSIN 00 0SN=F9G7LB. TPLOT. FORT, DISP=SHR 
//LKED, SYSPRINT UD DUMMY 
/ /GO . SYSIN OU * 

0.0(0.04 
0. 0(2.0 
0 o 0 ( 4 . 0 
0.0(100.0 
0.0,100,0 
“ISOoQsSO.O 
0 ,8 , 1 « 0 
0,8, I, 0 
0 c 0 « 1 • 0. 

0 o 0 ( 1 • 0 

1.0(1,25 

/# 

DO OSI\l=&aTRES,DXSP=< OLD, DELETE) 


<CF,F081) ,71X153015 « N0TIFY=F9G7L8 , 


//GO.FTllFOOl 
// 

//TPL0T3B JOB 
// TIME— 1 
/*ROUTE PRINT REMOTE2 
/*SETUP PLOTS, NOCOUE 

//TRES EXEC FORTH, LIB4= , IMSL, SINGLE* 
//FORT.SYSPRXNT OD DUMMY 

//FORT, SYSIN DO DSN-F9G7LB. TRES. FORT ,DISP=SHR 
// DD DSN=F9G7LB,M0DEL3,F0RT,DISP=SHR 

//LKED .SYSPRINT DD DUMMY 

//GO.FTOlFOOl UD DSN=F9G7LB» MQDEL38 .DATA , DISP-SHR 
//GO. SYSIN DO * 

1 o 0 , 1 , U , 0 . 005 

27,0.001 

f $ 

//GO.FTlOFOOl DO DSN=F9G7LB , DYN3B , DATA « D1 SP=SHR 
//GO.FTllFOOl UD DSN=&STRES,DISP=( .PASS) t 
// UNIT=DISK,SPACE=(TRK, (10,2) ) 

//LIST EXEC TSOLIST,NAME=*a&TRES' ,C0ND=EVEN 
//TPLOT EXEC FORTHP 
//FORT.SYSPRXNT DD DUMMY 

//FORT.SYSIN OU DSN=F9G7LB , TPLOT , FORT , DISP=SHR 
//LKED. SYSPRINT DD DUMMY 



151 


//GGoSYSIN DU * 

0,0,0,04 
0 q 0 f 5 4 0 

0,0,10o0 

0,0,500,0 
0,0,500,0 
*»500 , 0 , 0 , 0 
0 * 7 , 0 * 9 

0 « 8 « 0 » 9 

•2»oa«o 
-3 • 0 1 1 , 0 
X , 0 , 1 0 5 
/# 

//GO.FTUFOOl DD DSN=&£TRES , DISP= ( OLO, DELETE) 

//TPL0T4 DOB (CF,F081) ,7111530 15, N0TXFY=F9G7L8, 
// TII"IE=1 

/*ROUTE PRINT REM0TE2 

/#SETUP PLOTS, NOCODE 

//IRES EXEC FORTH, LIB4=» XMSL, SINGLE* 

//FORT.SYSPRINT DD DUMMY ■ 

//FORT ©SYSIN DO DSN=F9G7LB»TRES,F0RT ,DISP=SHR 
// DO DSN=F9G7LB,MOOEL4,FORT,DISP=SHR 

//LKED oSYSPRlNT DD DUMMY 

//GO , FTOIFOOI DD 0SN=F9G7LB , M00EL4 , OATA , DISP=SHR 
/ /GO«SYSIN DD * 

1 o 0 ,1 » 0 , 0 o 005 
800,0,001 

/ * 

//GOeFTlOFOOl DD DSN=F9G7LB, DYN4,DATA,DISP=SHR 
//GOoFTllFOOl DO DSN=&&TRES , DXSP= { , PASS ) , 

// UNXT=DISK,SPACE=(TRK, U0,2> ) 

’//LIST EXEC TSOLS ST » NAME= * &&TRES * , COND^EVEN 
//TPL'OT EXEC FORTHP 
//FORT«SYSPRINT DO DUMMY 

//FORT « SYS IN DO DSN=F9G7LB , TPLOT ,FORT , DI SP=SHR 
//LKED eSYSPRXNT DD DUMMY 
//GOoSYSIN DU * 

0 « 0 , 0 , 8 
0 « 9 « 1 » 0 
0,9, 1*0 
0 , 8 , 1,0 
0,85,1,1 

0 ♦ 8 , 1 o 0 

0 , 7 , 1 e 0 
0 # 8 , 1 e X 
0,85, Xal 
0 , 9 e 1 o 1 

0,85,1,1 

/* 

//GOoFTllFOOl DD DSN=£&TRES , DXSP= < OLD , DELETE ) 
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THIS PROGRAM COMPUTES THE TIME NESPONSE OF A MODEL 
UNDER THE CONTROL OF AN OPTIMAL CONTROL LAW « 

AS GENERATED 8Y THE DYNANIC PROGRAMMING METHOD, 

REAL X<5) ,DX(5} ?U<2> , Y<3) 

COMMON XMIN ( 2 ) ,XMAX(2) ,XXNC<2) ,XTAR<2> 

COMMON VOPT(15?15) ?UGPTC15»15?2) 

1 = 0,0 

CALL INIT(MDL) 

READ ( 10 ) IT? VOPT ?UOPT 

READ ( 10 ) XMIN, XMAX? XINC ? XTAR 

WRI TE ( 11 , 200 ) MDL 

READ ( 5 « * ) U ? EPS 

CALL STDST { U « X ) 

CALL CNTRL { X? U , VI ) 

WRITE ( 6 ? 30 0 ) VI 
WRITE ( 11 ♦ 201 ) T » X 
CALL OUTPUT ( X , U , Y ) 

WRITE (11,202) U ? Y 
5 READ(5?*,END=999) N,DT 
DO 10 1=1, N 
CALL CNTRL ( X ? U , V ) 

CALL XDOT (X?U,DX) 

UO 20 J = 1 » 5 
X( J)=X< J>+DX( J)#DT 
20 CONTINUE 

CALL OUTPUT ( X , U , Y ) 

T=T+DT 

WRITE (11,201) T ? X 
WRITE (11,202) U » Y 

D.IST=SQRT( (X(l)-XTAR(l) ) **2+ ( X ( 2 ) -XTAR ( 2 ) )**2) 

REL=(T»VI)/T 

I F ( DX ST , LT , EPS ) WRITE<6,301) T,KEL 
1F( DIST ,LT ,EPS ) EPS=0,0 
10 CONTINUE 
GO TO 5 
999 CONTINUE 
STOP 

200 FORMAT ( A4 } 

201 FORMAT ( El 0 , 3, 5E1A , 6 ) 

202 FORMAT (10X,5E14,6) 

300 FORMAT (IX, * PREDICTED TIME =««F8,4> 

301 FORMATdX, ’ACTUAL TIME =« ?F8,4/1X, ’RELATIVE ERROR s*«Ell*4> 
END 



XHAX(l)) GO TO 10 
XMIN(X) > /XINC { 1 ) + l o 0 


SUBROUTINE CNTRL ( X * U « V ) 

THIS ROUTINE INTERPOLATES CONTRuLS FROM THE 
ENTRIES OF THE OPTIMAL CONTROL LAW. 

. REAL X { 5 ) , U ( 2 ) , 

COMMON XMIN ( 2 ) ,XMAX{2> % XINC ( 2 ) , *TAR ( 2 ) 
COMMON VOPT ( 15 « 15 ) tUOPT (15,15,2) 

NX = 15 
XX1L=NX 
IX1H=NX 
FR1L=1.0 
FR1H=0.0 
IF ( X ( 1 ) *GE 

FRlHs(XU) 

IX1L-IF IX(FRIH) 

FR1H=FR1H»FL0AT(IX1L> 

IX1H=1X1L+1 
FR1L=1.0«FR1H 
10 CONTINUE 
IX2L-NX 
IX2H=NX 
FR2L=1 » 0 
FR2H=0.0 
IF ( X ( 2 ) oGE. 

FR2H~ { X ( 2 ) < 

IX2L=IFIX(FR2H) 

FR2H=FR2H~FL0AT(XX2L) 

XX2H=IX2L-s-l 
FR2L=X.0-FR2H 
20 CONTINUE 

V=FR1L*FR2L*V0PT(IX1L,IX2L> 
U{1)=FR1L*FR2L*U0PT(IX1L,IX2L»1) 
U(2)=FR1L*FR2L*U0PT{ IX1L»IX2L»2) 
V=V+FR1H*FR2L*V0PTUX1H,XX2L> 
U(1)=U‘(1)+F.R1H*FR2L*U0PT( IXlH«Ix2Ltl) 
U(2)=U(2)+FR1H*FR2L*U0PT( IX1H«IX2L*2) 
V=V+FR1L*FR2H*V0PT (XxXL, IX2H) 
U(11=U(1)+FR1L*FR2H*U0PTUX1L«IX2H»1) 
U(2)=U(2)*s-FRlL*FR2H*UQPT ( I XXL. Ix2H» 2) 
V=V+FR1H*FR2H*V0PT(IX1H»IX2H) 

U(1>=U( 1>+FR1H*FR2H*U0PT< IxIHi I x2H,l) 
U(2)=U(2)4-FR1H*FR2H*U0PT( XXlH, Ix2H,2) 
RETURN 
END 
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, XMAX ( 2 ) ) GO TO 20 
•XMIN ( 2 ) )/XINC(2)4-1 o 0 


SUBROUTINE OUTPUT < X , U , Y } 

THIS ROUTINE o COMMON TO ALL MODtLSi 
EVALUATES OUTPUT VARIABLES. 

REAL X<2) «U(2) *Y(3> 

REAL C ( 3 , 2 ) ? D ( 3 » 2 ) « E ( 3 ) , WK ( 3 ) 

DATA C/-0. 61059, -0. 20154, -0.582*9,-0. 10759. -0. 45613 « 0,46872/ 
DATA D/0, 50292, 0.20 423 ,0.1 8877, u. 17669 , 0,14724 ,-0.92545/ 

DATA E/1,03837, 1.30796, 1.85021/ 

CALL VMULFF (C,X»3,2»1,3»2»Y,3,ILR) 

CALL VMULFF (D,U, 3,2, 1,3,2, WK,3,iER) 

DO 10 1=1,3 
Y (I)=Y(I)*WK( I)+E{I) 

10 CONTINUE 
RETURN 
ENO 



APPENDIX E 


TIME RESPONSE PLOT PROGRAM 

This appendix includes the time response plot program, designed to plot 
ten graphs: states, controls, and outputs versus time. The JCL in appendix 

D links this program to the general system simulator, and the JCL in appen- 
dix A links it to the DYNGEN simulator. In either case, this program’s in- 
put consists of a data set of time response values, and user specified axis 
limits . 
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THIS PROGRAM PLOTS THE TIME RESPONSE 
DATA GENERATED BY TRES, INCLUDING 
STATES, controls. AND OUTPUTS. 


REAL TS ( 1003 ) ,XS(1003.5) . US ( 1003 , 2 ) i YS ( 1003 « 3 ) 
INTEGER ITXT(2.10) 

DATA ITXT/*X(I) *«'=NC * . 9 X ( 2 > 9 * 9 =NF 9 . 

S*X<3) 9 , 9 =P4 * , 9 X(4) * . *=P7 « , 9 X( 3 ) 9 , 9 =U4 * « 


& 9 Ud) 9 , 9 =WFB 9 . 9 U(2) 9 . 9 =A8 ». 

&*Y(1) 9 i »=:ZC 9 . 9 Y(2> * . 9 =ZF 9 « 9 Y ( 3 ) * * 9 =T4 


V 


TLEN=Q s 0 
XLEN=5 « 0 

READ (11.200) MDL 
DO 2Q NS=1 *1001 

READ(11.*.END=21) TS ( NS ) « ( XS l NS ♦ J ) tjslt 5) 

READ<11***END=2X> ( US ( NS , J ) .J=1.2) . ( YS ( NS. J ) * J=1 ♦ 3 ) 

20 CONTINUE 

21 CONTINUE 

NS=NS-1 ■ s 

-CALL PLOTS (711153015) 

CALL PLOT (0.0.0s5.3) 

CALL PLOT (0«0*10«5.2) 

CALL PL0Tda75.l85.-3) 

READ15.*) TMIN.TMAX 
TS(NS+X)=TMIN 
TS ( NS+2 ) = ( THAX-TMIN ) /TLEN 
DO 30 J=l*5 

- CALL SYMBOL(-Oe75(O e 5.0*l4. 29H CONTROLLER RESPONSE OF MODEL 
&90o0.29> 

CALL SYMBOL (999,0. 999,0. 0.14 .MDL .90,0 .4) 

CALL SYMBOL (-0«5.0*5.Q»14« 23HF0K INITIAL X=U. 0*1,0). 

&9Q o 0 .23 ) 

READ < 5 » # ) XMIN.XMAX 

XS(NS+1*J)=XHAX 

XS(NS+2*d)=(XMIN-XMAX)/XLEN 

CALL AXXSC0,Q*0,0»ITXT(l*d} « -8 . XLEN » 0 . 0 * 

#XS(NS + 1, J) « XS ( NS+2 . J ) ) 

CALL AXIS(XLEN*0.0»15HTIHE IN SlCONDS « - 15 « TLEN . 90 , 0 • 
8TS(NS+l)«TS{NS+2) ) 

CALL LXNE(XSd.J) ,TS*NS*X* 0*0) 

CALL PL0T(6o75.-lo0.3) 

CALL PL0T(6,75. 9,0.2) 

CALL PLOT (085.0sC.~3) 

30 CONTINUE 

CALL°SYMb6l(-0,75,0,5,0,X4.'29HCuNTROLLER RESPONSE OF MODEL 
&9Q e 0 . 29 ) 

CALL SYMBOL ( 999 o 0.999, 0.0 . 14.MDL .90 a 0.4) 

CALL SYMBOL('*0o5. 0,5.0, 14. 23HF0K INITIAL X=(X,Q*l.Q>* 

&90 » 0 , 23 )• 

READ ( 5 . # ) UMIN » UMAX 
US(NS-d.U)=UMAX 


. 
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■ tt « XLEN » 0 o 0 * 


US { NS4-2 « J ) = ( UMIN°UMAX ) / XLEN 
CALL AXXS(0,O«0c 0»1TXT(1« J-J-5) t' 

8US(NS + 1<) J) *US<NS + 2» U) ) 

P, f 9 ♦ 15HTXME IN StCONDS»-15«TLEN,90»0, 

« 1 S ( NS+1 ) 9 I S ( NS + 2 ) ) 

CALL LIN£(US(1, J) « TS«NS« 1 « 0 « 0 ) 

CALL PLOT(6.75»-1.0t3) 

CALL PLOT (6s75t9«,0»2) 

CALL PLOT <8 P 5«0«0»“3) 

40 CONTINUE 
DO 50 d=l,3 

|YFlBOL('*0 c 75«0 9 5»0 o l4 s a9HCuNTROLLER RESPONSE OF MODEL 
« V U « U vt?) 

CALL SYMBOL ( 999 , 0 * 999 ♦ 0 * 0 » 14 1 MOL ♦90, 0*4) 

CALL SYMBOL{°0«5<0*5»0c14«23HFOK INITIAL X= { 1. G , 1 » 0 > , 

&y u « u « 2 o ) 

READ (5s*) YMIN « YMAX 

YS(NS+1,U>=YHAX 

YS(NS+2« J)s<YHIN-YHAX)/XLEN 

„S§bbe A ?^! 0 e2^ 0 A°l ITXT{:L «' J+7 >’‘' o ’XCEN»0 e Q» 

&^S(NS+1?U) 9YS(NS+2»J) ) 

«?cVhc^?? S 42V5S*??9 ,15HTIME: IN SECONDS, -15, TLEN, 90,0, 
tfTS ( NS+1 ) » TS ( NS +2 ) ) 

CALL LINE<YS<l.d) «TSiNSil t Q»0> 

CALL PLOT (6»75,-l«0,3) 

CALL PLOT (6.75?9c0t2) 

CALL PLOT (8.5»0»0»**3) 

50 CONTINUE 

CALL PL0T{Q e Q,0.0,999) 

RETURN 

200 FORMAT (A4) 

END 



APPENDIX F 


DYNAMIC SUCCESSIVE APPROXIMATION PROGRAM 

This appendix includes the JCL used to generate control laws, and the 
dynamic successive approximation program. The program is designed to be 
compiled with a model program consisting of subroutines named INIT, XDOT, 
OUTPUT, and COMPLY. The program’s input includes state and control space 
limits and quantization, target point, time increment, and output constraint 
values. A data set containing the control law (initialized to the target’s 
steady state controls) and associated cost estimates is read by the program. 
A user specified number of iterations is performed, and the new control law 
is rewritten to the input data set. Subsequent runs then perform more 
iterations on this updated data set. 
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X5t> 

//DYNPRG3A JOB ( CF , F081 ), 7X1153015 « NOT IF Y=F9G7LB , 

// TIME=4 

/♦ROUTE PRINT REM0TL2 

//DYNPRG EXEC F0RTH«LIB4= f XMSLe SINGLE* 

//FORToSYSPRINT DO DUMMY 

//FORT o SYSIN DU DSN=F9G7LB • DYNPRG* FORT , DXSP=SHR 
// DD DSN=F9G7LB*M0DEL3,F0RT,DISP=SHR 
//LKEDoSYSPRINT DD DUMMY 

//GO.FTOIFOQI DD DSN=F9G7LB , M0DEL3A . DATA , DISP=SHR 
//GO'« SYS IN DD * 

0*948434*0, 928751 
0 , 01 * 0,01 
0,74*0,74 
1 , 0 , 1,0 
0 , 02 , 0,02 

1,15,1,105,1,08 

5. 250. 0. 01.0.01 
T,F»F 

/♦ 

//GOoFTiOFO 0 1 DD DSN=F9G7LB , DYN3A , DATA , DISP=OLD 
// 

//DYNPRG3B JOB (CFtFOQl ) , 71H53015 , N0TIFY=F9G7LB , 

// TXME=4 

/♦ROUTE PRINT REM0TE2 

//DYNPRG EXEC FORTH » LI B4= * IMSL , SINGLE * 

//FORToSYSPRINT DD DUMMY 

//FORT o SYSIN DD OSN=F 9G7LB . OYNPRG ,F0RT , DISP=SHR 
// DD OSN=F9G7LB» MODEL 3 , FORT , DI SP~SHK 
//LKED,SYSPRXNT OD DUMMY 

//GOoFTOlFOOl DD DSN=F9G7LB, M0DEL3B , DATA , DISP=SHR 
//GOoSYSlN DU * 

0*948434,0,928751 

0 , 01 , 0,01 

0,74,0,74 

1 . 0 . 1 . 0 
0 o 02 , 0,02 

1,15,1,105, 1,08 

5.250.0. 01.0.01 
T , F , F 

/# 

//GOoFTlOFOOl DD DSN=F9G7LB , DYN3B « DATA , D1SP=0LD 
// 

//DYNPRG4 JOB ( CF , F081 ), 7111530 15 , N0TIFY=F9G7LB , 

// TIME=9 

/♦ROUTE PRINT REMOTE2 

//LIST EXEC T SOL I ST , NAME= * F9G7LB , M0DEL4 , DATA* 

//OYNPRG EXEC FORTH, L1B4=® I MSL« SINGLE* 

//FORToSYSPRINT DD DUMMY 

//FORT, SYSIN DU DSN=F9G7LB, DYNPRG, FORT, DISP=SHR 
// DD DSN=F9G7LB»M0DLL4oF0RT,DISP=SHR 
//LKEDoSYSPRINT OD DUMMY 

//GOoFTOlFOOl DD DSN=F 9G7LB , MQDEL4 ,QAT A , DISP=SHR 
//GOoSYSlN DD * 

0,94 8434,0® 928751 

0.0i,0o0i 

0,74,0,87 

1.0. 1*13 

0 , 02 , 0,02 

1,15,1,105,1,08 


1 , 250 , 0 , 01 , 0,01 

T , F , F 

/t- 

//GOoFTlOFOOl DD DSN=F9G7LB , DYN4 , DATA, DISP=0LD 
// 
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THE FOLLOWING VARIABLES ARE USED IN THIS PROGRAMS 


•MDL 

NX 

IXTAR 

IWStZWE 

NNX 

XTAR 

XINC 

IX 

XS 

XMIN 

XMAX 

UMIN 

UMAX 

UINC 

NU1 » NU2 

IU1»IU2 

YMAX 

NT 

NTS 

DT 

DTI 

IT 

VOPT 

UOPT 

FRIT 

PRITS 

RUN 

XXI, 1X2 

NCHNG 

XIN 

XX NS 

XOPT 

ITS 

NCHNGS 


MODEL NUMBER (ALPHANUMERIC) 

STATE SPACE DIMENSION 
TARGET POINT INDEX 

STATE SPACE WINDOW START AND END INDICES 

NUMBER OF STATE SPACE POINTS EXCEPTING THE TARGET 

TARGET POINT VALUES 

STATE SPACE INCREMENT 

STATE SPACE INDEX “ 

STATE SPACE VALUES 
STATE SPACE MINIMA 
STATE SPACE MAXIMA 
CONTROL SPACE MINIMA 
CONTROL SPACE MAXIMA 
CONTROL SPACE INCREMENT 
CONTROL SPACE DIMEi\SlONS 
CONTROL SPACE INDICES 
OUTPUT CONSTRAINT MAXIMA 
MAXIMUM ITERATIONS 
MAXIM! UM SUBITERATIONS 
TIME INCREMENT 
INTEGRATION TIME INCREMENT 
ITERATION NUMBER 
COST ESTIMATES 
OPTIMAL CONTROLS 

ITERATION PRINT CONTROL (LOGICAL) 

SUBITERATION PRINT CONTROL (LOGICAL) 

RUN CONTROL (LOGICAL) 

STATE SPACE INDICES 
UPDATE COUNTER 

ALLOWABLE VALUE XNUICATOR (LOGICAL) 

ARRAY OF SAVED XIN*S 
ARRAY OF SAVED RESPONSES 
SUBXTERATION number 
UPDATE COUNTER 


REAL UMIN ( 2 ) ,UMAX(2) ,UINC(2> »US(21,2) 

REAL X ( 2 ) « XN ( 5 ) «U<2) 

LOGICAL XIN, PRIT, PRITS* RUN, XINS(I5» 15) 

COMMON DT,DTI«XXl*IX2«NX,XIN*XXfAR c ^ v 

COMMON XMI-N ( 2 ) , XMAX { 2 ) »XXNC(2) « xTAR ( 2 ) , XS U5 « 2 ) ,YMAX(3> 
COMMON VOPT (15,15) , UOPT ( 15 * 15 , 2 ) , XOPT ( 15 « 15 » 2 ) 

CALL INXT(MDL) 

NX=15 

IXTAR=(NX+l)/2 
!WS=IXTAR-5 
IWE=I XTAR+5 
NNX=NX$*2«i 
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define: the elements and dxmensiuns 

OF THE STATE SPACE 

READ<5t*) XTAR 
WRITE ( 6 9 220 ) XTAR 
READ ( 5 « ^ } XINC 
WRITE ( 6 ? 221 > XINC 
DO 10 IX=1,NX 

XS ( XX a ) =XTAR (1 H-XINC < 1 > *FUOAT< AX-IXTAR ) 
XS( IX,2)=XTAR(2)+XINC<2)*FLQATUX-IXTAR> 
10 CONTINUE 

XKIN<l)sXSUti> 

XNIN(2)=XS(i92) 

XMAX { 1 ) sXS ( NX « 1 ) 

XMAX ( 2 ) ~XS ( NX ♦ 2 } 

DEFINE THE ELEMENTS AND DIMENSIONS 
OF THE CONTROL SPACE 

READ { 5 » * ) UMXN 
WRITE { 6 v 222 ) UMIN 
■ READ { 5 » * ) UMAX 
WRITE ( 6 9 223) UMAX 
READ ( 5 9 * ) UINC 
WRITE { 6 9 224 ) UINC 

NU1~IF X X { < UMAX U) “UMXN U ) > /UINC U ) +0 *5 > +1 
DC 20 IU1 = 1 1 NU1 

US( IU1»1>=UMIN< D+UINC U) *FLOAT UU1-1) 

20 CONTINUE 

NU2=IFIX( <UMAX(2)«UMIN<2> > /UINC < 2 > +0 ,5) +1 
DO 2b IU2=1 « NU2 

US( XU2»2>=UMIN{2)+UINC(2)#FL0ATUU2-1) 

25 CONTINUE 


INPUT OTHER CONTROL PARAMETERSAnD 
INITIAL COST ESTIMATES 

READ ( 5 » « ) YMAX 
READ ( 69 *) NT9NTS»DT,DTI 
READ< 10 ) XT 9 VOPT 9 UOPT 
WRXTE<6,225) NT 9 NTS 9 DT 9 DTI 9 IT 
REWIND 10 

READ ( 5 « * } PRITtPRITStRUN 
VOPT( IXTAR 9 IXTAR)=0.0 
IF ( ® NOT 0 RUN ) GO TO 90 


SCAN THE STATE AND CONTROL SPACES* 
EVALUATING THE COST 'FOR EACH SE ! OF 
CONTROLS « REPLACE THE CONTROLS IF 


THE COST IS LESS THAN THE 
INITIAL ESTIMATE* 


DO 70 vJT=i » NT 

XT-IT-s-1 

1X1 = 1 XT AR 

IX2=XXTAR 

NCHNG=0 

DO 60 IX=1«NNX 



CALL SPIRAL 
Xd)sXS(IXltl) 

X(2)=XS(IX2 9 2) 

XINSC 1X1,1X25=, FALSE, 

DO 50 IU1=1»NU1 
DO 50 IU2=1«'NU2 
U(l)=USUUl # l) 

U { 2)=US ( IU2 ( 2 } 

CALL NEXTX ( X ? U ? XN ) 

IF { ,N0T oXIN>GO TO 50 
VTST-V ( XN) 

IF ( VTST aGE<>V0PT(IXl«IX2) ) GO TO 50 

NCHNG=NCHNG+1 

V0PT(IX1»IX2)=VTST 

U0PT(IX1 5 IX2»1)=U(1> 

U0PT(IX1«IX2»2)=U(2) 

XOPT ( XXI « XX2«1)=XN(1) 
X0PT(XX1,IX2»2)=XN(2) 
XINS(2X1,IX2)=XIN 
50 continue; 

60 CONTINUE 

IF(PRXT) WRITE 1 6 « 210 ) IT * NCHNG 
IF (NCHNG«EQo 0 ) IT = IT-1 
IF ( NCHNGe EGU 0 > GO TO 80 

SCAN THE STATE SPACE ONLY « 
REEVALUATING THE COST FOR 
A FIXED CONTROL LAW, 

ITERATE UNTIL THE TRUE COST 
FOR THIS LAW IS FOUNO, 

DO 40 ITS=1»NTS 

NCKNGS=0 

IX1=IXTAR 

XX2=IXTAR 

DO 30 IX=1,NNX 

CALL SPIRAL 

IF( ,NOT,XINS( XXI? 1X2) ) GO TO 30 

XN.(l>-sXOPTUXltIX2tl> 

XN(2)=X0PT(XX1*IX2»2> 

VTST=V(XN> 

■IF<VTST*GE,V0PT<IXX,.IX2) > GO TO 30 
NCHNGS=NCHNGS+1 
VOPT ( IX1»IX2)-VTST 
30 CONTINUE 

IF(PRITS) WR ITE ( 6 » 211 ) ITS? NCHNGS 
IF(NCHNGSoEQ«0) GO TO 45 
40 CONTINUE 
45 CONTINUE 
70 CONTINUE 
8U CONTINUE 

RECORD ANO REPORT THE RESULTS 

WRITE 1 10 I I T « VOPT o UOPT 
WRITE (10) XHIN*XMAX«XINC*XTAR 
90 CONTINUE 

WRITE(6«200» 
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100 


110 


WRITE (6,201) {XSlIXhlh XX1=1WS»IWE) 

DO 100 IX=IWS,IWE 

WRIT£*6*202) XS ( 1X2,2 ) , ( VOPT ( IXl , 1X2 ) , 1X1=1 WS , I WE ) 
CONTINUE 

IF( YMAX(1)*YMAX(2)*YMAX(3) .EQ.O.Q) WRITE (6, 20 7) 

IF{ YMAX(1)*YMAX(2)*YMAX(3) .NE.0.0) WRITE (6,208) 

WRITE ( 6 , 200 ) 

WRITE (6, 201) < XS( 1X1,1) , IX1=1WS,IWE) 

DO 110 IX=IWS,IWE 
IX2=NX«- I X-S-l 

XS{ 1X2*2) , (UOPTC IXl « 1X2,1) , XX1=XWS, 1WE> 
(UUPT(IX1,IX2,2) , Xxl=IWS» IWE ) 


120 

200 

201 

202 

203 

204 

205 

206 

207 

208 
210 
211 
212 
213 
220 
221 
222 

223 

224 

225 


WRITE (6*203) 
WRITE ( 6 * 204 ) 
CONTINUE 
WRITE (6, 206) 


MDL 


WRITE ( 6 » 207 ) 
WRITE (6,208) 


IF ( YMAX ( 1 ) *YMAX ( 2 ) *YMAX { 3 ) , EQ.O « 0 )’ 
XF(YMAX(1)*YMAX(2)*YMAX(3) .NE,0«0) 

WRITE ( 6 , 200 ) 

WRITE { 6 * 201 ) {XS (1X1,1) , 1X1 = 1 WS, IWE ) 

DO 120 I X=I WS , IWE 
IX2=NX“' IX + 1 

WRITE (6, 212) XS(IX2,2) , (XOPT( IXl ,1X2,1) , 1X1=1 WS, IWE) 
WRITE (6,213) (XOPT( 1X1, 1X2, 2) , IXl=I WS, IWE > 

CONTINUE 

STOP 

FORMAT ( ’ !♦///////////////) 

FORMAT (1X«T21*11F8.4///) 

FORMAT (/8X,F8 e 4,T21,llF8o4/) 

FORMAT (/8X,F8 e 4,T20,llF8»2) 

FORMAT (1X,T20«11F8.2) 

FORMAT ( ///IX, T20, » FIGURE 
FORMAT ( ///IX, T20 , ’FIGURE 
& * OPT IMAL CONTROL LAW*) 

FORMAT ( ♦ + » * T50 , 'UNCONSTRAINED* ) 

FORMAT ( * + • , T50 , ! CONSTRAINED* ) 

FORMAT ( IX, T5 0, 5 ITERATION = * ,X4, 

FGRMATUX, ’PREITERATION = «, 
FQRMAT(/8X,F8.4,T21,11F8«4) 

FORMAT ( 1X,T21,11F8,4) 

FORMAT (IX, 5 XT AR= { ' ,F8.4, 

FORMAT (IX, * XINC=( 

FORMAT (IX, *UMIl\I= ( 

FORMAT ( IX , * UMAX= ( 

FORMAT (IX, *UINC=( 


■A 

>6 


MODEL 

MODEL 


♦ , A4,T65, ’COST* ) 
* , A4,T65, 


14 < 


, NUMBER OF UPDATES =*,I5) 
NUMBER OF UPDATES =*?I5) 


,F8#4, 

» F8 . 4 , 

« F8 . 4 , 

? F0 q 4 t _ 

FORMAttlxi *NT=*,i3/lX, * NTS= * , I 4/ IX , 5 DT= * ,E11,4/ 
&1X, *OTI=* ,Ellc4/lX, , *IT=* ,13/) 

END 


»F8 e 4 , * > 
, F8 » 4 , * ) 
, FS » 4 » * > 
, F8 , 4 , * ) 
, F8 » 4 , * ) 


10 

20 

30 

40 


SUBROUTINE SPIRAL 

THIS ROUTINE UPDATES THE 
STATE SPACE INDICES , SPIRALLING 
OUTWARD FORM THE TARGET. 

COMMON 0T«0TI*IXl,lX2*NX«XINtIXlAR 
COMMON. XMXN ( 2 ) ,XMAX<2> »XXNC(2) , xTAR ( 2) ,XS( 15,2) ,YMAX(3) 
COMMON VOPT (15,15) , UOPT ( ._ I . 1 1 
IF{ XXI* SX2.GE,2» X XTAR, AND ,1X1 
IF ( 1X1+ 1X2. GE* 2*1 XT AR. AND. 1X1 
IF ( IX1 + IX2.LT »2*IXTAR» AND. 1X1*6) 

XF( IXl+IX2,LT,2*IXTAR.ANDt IXl*Lt. 

STOP 16 
IX!=IXl+i 
RETURN 
IX2=iX2«l 
RETURN 
IX1=IX1-1 
RETURN 
1 X2=I X2+1 
RETURN 
END 


,15 

,2) 

, XOPT 

( 15 < 

>15< 

.2) 

1X1 

« LI 

.1X2) 

GO 

TO 

10 

1X1 

• Gt 

. 1X2) 

GO 

TO 

20 

1X1 

. G ) 

,1X2) 

GO 

TO 

30 

XXI 

«Le. 

.1X2) 

GO 

TO 

40 



163 


20 


30 


40 


REAL FUNCTION V(XN> 

THIS ROUTINE EVALUATES THE COST 
FUNCTION AT ANYPOINT IN THE 
STATE SPACE* BY INTERPOLATING 
FROM NEIGHBORING GRID VALUES* 

REAL XN(5) 

COMMON l ~OT«OTI*IXltXX2*NX«XIN*XX»AR 

COMMON XMIN<2> *XMAX<2> *XINC(2> *XTAR(2 *XS(15*2>*YMAX(3> 
COMMON VOPT ( 15*15) «UOPT (15*15*2) * XOPT ( 15 * 15 * 2 ) 
FR1H=(XN(1> -XMIN(l) ) /XINC ( 1 ) +1 * U 
IX1L-XF IX (FR1H ) 

FR1H=FK1H~FLQAT< IX1L) 

IXIH-IXIL+I 
FR1L=1 « 0-FR1H 

FR2H= (XN(2)“XMXN(2) ) /XINC { 2 ) +1 , U 
IX2L=IFIX(FR2H) t 

FR2H=FR2H“FL0AT(IX2L> 

IX2H=XX2L-J-1 
FR2L-1 * 0**FR2H 

V=QT*FRlLsFR2L*V0PT ( IX1L.IX2L) 

IF ( IXiHcGT.NX) GO TO 20 4 

V=V+FR3LH*FR2L*VOPT( IXIHt IX2L) 

CONTINUE _ , 

IF ( I X2H 0 GT « NX ) GO TO 30 
V=V+FK1L*FR2H*V0PT(IX1L*IX2H) 

XFaXlH»GT.NX e 0R.IX2H.GTeNX) GO TO 40 
V=V+FRlH*FR2H*V0PT(IXlHiIX2H) 

CONTINUE 

RETURN 

ENO 


SUBROUTINE NEXTX ( X » U * XN ) 

THIS ROUTINE PERFORMS A VARIABLE 
STEP EULER INTEGRATION, AND 

TESTS THE NEW STATE AGAINST ALL CONSTRAINTS, 

REAL X(2> *U<2) *XN(5) *DX<5) *Y(3> 

LOGICAL XIN 

COMMON DT,DTI*IX1*XX2*NX*XIN«IX1AR 

COMMON XMIN ( 2 ) ,XMAX(2) * XINC (2) , XTAR ( 2 ) * XS( 15 * 2 ) , YMAX (3) 
COMMON VOPT (15, 15) * UOPT <15 « 15 , 2 ) * XOPT ( 15 * 15 , 2 ) 

DO 10 J=l*2 
XN< J)=X( J> 

10 CONTINUE 

CALL CUMPLT ( XN * U } 

NI=XFIX(DT/DTI+G,5> 

DO 20 1=1, NT 
CALL XOQT ( XN * U * DX ) 

DO 20 J=l*5 

XN( J)=XN(J)+DTI*DX( J) 

20 CONTINUE 
XIN=»TRUE9 

IF ( XN ( 1 ) oLToXMINU) ) 

IF (XN(2) ,LT e XMIN ( 2 ) ) 

IF ( XN ( 1 ) *GT »XMAX ( 1 ) ) 

IF ( XN ( 2 ) <i GT , XMAX ( 2 ) ) 

CALL OUTPUT ( X , U * Y ) 

IF ( Y ( 1 ) ,GT* YMAX(l) ♦ AND » YMAX ( 1 ) .OT.OoQ) 

IF ( Y ( 2 ) c GT » YMAX ( 2 ) * AND, YMAX (2) ,bT,0„0> 

IF ( Y ( 3 ) 0 GT, YMAX (3) .AND, YMAX (3) ,bT,0.0) 

RETURN 
END 


XIN=, FALSE, 
XIN=. FALSE, 
XIN=, FALSE, 
XIN=, FALSE, 


XIN=. FALSE, 
XIN=. FALSE, 
XIN=, FALSE, 
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SUBROUTINE OUTPUT (X»U,Y) 

THIS ROUTINE i COMMON TO ALL MODELS , 

EVALUATES OUTPUTS FOR*CONSTRAlN \ COMPARISON 


X<2> »U(2) ,Y(3> 

C<3,2> cD(o»2) ,E(3) tWK(3) 


C/-0 « 610 59 , -Q, 20 15A, -0*58229, ~0„ 10759 * - 0 e 45813 « 0,46872/ 
D/0 s 50292, 0. 20423 « 0 , 16877, U , 17689 ♦ 0 , 14724 ,*» O . 9254 5/ 

tr S 1 A70 - 2’? -I ~Z ~7 Q S' 1 flCAOl / 


W 1 ' W t W •* >- 1 V W •— v * t » » * I 

DATA E/1,03637, 1,30796,1, 85021/ 

CALL VMULFF(C,X»3t2»l«3,2»Y»3.»ItR> 
CALL VMULFF<D,U*3*2,1*3»2«WK»3, iER-> 
DO 10 1=1,3 
Y (I)=Y(X)+WK(I)+E(r> 

CONTINUE 

RETURN 

END 




